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Rapid Creep in Structures 


N. J. HOFFT 
Polytechnic Institute of Brooklyn 


SUMMARY 
A survey of the effects of rapid creep upon the stress distribu 
tion in structures and of the various modes of creep failure is pre 
sented. It is proposed that rapid creep should be tolerated when 
aerodynamic heating causes high temperatures for only a short 
time in structural elements of supersonic aircraft 


THE CREEP PHENOMENON 


— OF THE ENGINEERING aspects of the 
creep phenomenon can be suitably begun by 
citing two papers published by Andrade just before 
World War I.' He observed that a weight suspended 
from the lower end of a vertical bar caused not only 
instantaneous, but also continued, elongations that 
increased with time. The tests were carried out with 
lead, tin, copper, iron, mercury, brass, and two other 
alloys at various temperatures. The history of the 
elongations following the application of the load in a 
typical test is shown in Fig. 1. 

Immediately after the load is applied, elastic and 
frequently also plastic deformations take place prac- 
tically instantaneously. As time goes on, the rate of 
elongation decreases rapidly during the so-called pri- 
mary or transient phase of creep. When the constant 
rate of elongation represented by the straight line is 
reached, the phenomenon enters its secondary phase. 
The slope of the straight line is the steady, or mini- 
mum, rate of creep whose values have been determined 
by many investigators for different metals at various 
temperatures and stress levels. Finally, in the ter 
tiary phase of creep, elongations develop at an increas 
ing rate until the bar fractures. 
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Metal physicists have proposed a number of dif- 
ferent mechanisms to explain the deformations shown 
in Fig. 1. The instantaneous elastic elongation is a 
consequence of changes in the lattice distance in the 
large number of crystals of which metals are composed 
and is recoverable upon removal of the load. When 
the tensile stress caused by the load exceeds the elastic 
limit of the material, plastic deformations also take 
place almost instantaneously as the load is applied. 
The mechanism of the plastic deformations is the slip 
process in which the material slides along favorably 
oriented planes like a pack of cards. 

Orowan* considers the primary phase of creep also 
as a consequence of this slip process. On the other 
hand, he attributes the secondary creep to continuous 
deformations taking place in the noncrystalline layers 
that fill in the small gaps between the crystals. These 
deformations are of a viscous type, resembling the 
motion of thick oil; as they take place, the adjacent 
crystals slide and roll. 

Both the primary and the secondary parts of the 
creep deformations are nonrecoverable. The deforma 
tions remain when the load is removed. The instan- 
taneous plastic deformation at ¢ = 0 cannot be recovered 
either; upon removal of the load only the elastic part 
of the instantaneous deformation is regained. This is 
usually small as compared to the permanent deforma 
tions when the stress and the temperature are high 
and when the load is applied for a considerable length 
of time. 

Andrade has shown that the tertiary phase of creep 
is a consequence of the test setup. It exists largely 
because of the particular properties of deformation 
geometry and load application. This tertiary phase, 
and the consequent fracture, will be discussed later 
at some length. At this time it need only be mentioned 
that under different conditions of loading the tertiary 
phase of creep can consist of decelerated motion. This 
is illustrated in Fig. 2, which shows the results of a 
bending creep test carried out recently at the Poly 


technic Institute of Brooklyn 
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Results of creep tests are usually presented in dia 
grams in which the applied stress is plotted on double 
logarithmic paper against the steady (secondary) creep 
rate and against the time necessary for fracture to 
With most materials the lines connecting the 
Thus the em- 


oecur.*~° 
test points are approximately straight. 
pirical relationship can be deduced 


€ = de/at = (e/r)" (1) 


where € is the steady creep rate, / is time, o is the tensile 
stress, \ is a constant, and m is a number equal to the 
slope of the straight line plot of stress versus creep 
rate. This type of plot as well as the power law were 


probably first used by Norton.‘ 


A few typical values of \ and m are given below. 
They were determined from data presented in a paper 
by Dorn and Tietz.’ The values are not very accurate 
because the basic information was taken from small- 
scale diagrams and because the linear approximation 
was not always satisfactory in the entire range of the 


stress. 


CREEP IN THE STRUCTURAL PARTS OF MISSILES 


The engineering theory of creep was originally de- 
veloped in order to analyze the deformations of rotating 
discs of gas turbines and the loosening of bolts in ma- 
chinery operating at high temperatures, as well as other 
problems of power engineering.*~'* These phenomena 
take place slowly; for this reason the theory had to 
predict deformations after thousands or ten-thousands 
of hours. 


The situation is entirely different with guided mis- 
siles whose duration of flight is measured in seconds or 
at most in minutes. Even if new power plants which 
permit longer flights at very high speed become avail- 
able, the duration of flight time will remain short. It 
should be remembered only that a missile flying at a 
Mach Number of 20 can reach any point on the surface 
of the earth in an hour. 


TABLE | 
Values of Exponent » in Creep Law from Tests by Dorn and 
Tietz 

Material Temperature 90°F, 212°F 300°F 100°F 
3S-H12 216 115 34 12.3 
35-H18 120 18 18.3 7 
528-H32 16.6 9.7 6.2 
52S-H38 38.7 13 6.0 5.2 
61S-T6 112 90 38.2 10 

248-T3 30 14 8.5 

TABLE 2 
Values of Constant A in Creep Law from Tests by Dorn and 
Tietz 

Material Temperature 90°F 212°F. 300°F $00°F 
3S-H12 19,600 17,600 15,800 12,500 
3S-HI1S 28,200 24,800 23,700 15,700 
528-H32 $3,000 38,000 22 000 
52S-H38 54,000 54,000 9 54,000 25,000 
61S-T6 19,000 42,000 42,000 46,000 
24S-T3 81,000 100,000 71,000 
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As missiles are almost always expendable, their 
design must be developed from the consideration that 
their structural parts need not remain strong and stiff 
beyond a few minutes of operation. Thus creep rates 
of several per cent per hour can be tolerated. Thy 
criteria are that the deformations developing during 
the design lifetime of the missile must not be large 
enough to interfere with the proper functioning of the 
moving parts of the missile, that the aerodynamic 
properties of parts exposed to the air stream must not 
be altered unduly, and that the structure must not 
fracture under the loads imposed in flight. 

The treatment of creep in missiles differs therefore 
considerahly from that developed for the purposes of 
mechanical engineering. Substantially larger defor- 
mations can be permitted in missiles, with the conse- 
quence that the initial purely elastic deformation is 
often reduced to a rather insignificant portion of the 
total deformation. On the other hand, at the same 
temperature higher stresses are allowed, and thus the 
instantaneous plastic deformation and the transient 
phase of creep gain in importance. (For structural 
design considerations in supersonic aircraft, see refer- 


ence 13.) 


FRACTURE IN TENSION 


When a rod is subjected to tension, it elongates and 
at the same time its cross section contracts. This 
phenomenon was discovered by Poisson, after whom 
the ratio of linear contraction to elongation is named." 
It has been customary to assume that plastic and vis- 
cous deformations take place without any change in 
the volume of the material. Consequently, Poisson's 
ratio is 1/2 when the deformations are small. 

Creep tests are usually carried out by suspending a 
weight from a vertical rod and leaving it there until 
the rod fractures. As long as the elongations are 
small, the cross-sectional area changes but little; and 
the tensile stress in the rod, which is the load divided 
by the area, remains sensibly constant. But when the 
area is decreased noticeably, the tensile stress increases 
appreciably. Such an increase in stress is accom 
panied by an increase in the strain rate; this latter 
can be considerable when v is a large number. 

The result is that the stretching of the rod is accel 
erated until the cross-sectional area is decreased to such 
an extent that it can no longer support the applied 
This is the explanation of 

The correctness of the 
He suspended 


load and the rod fractures. 
the tertiary phase of creep. 
explanation was proved by Andrade. ' 
a shaped weight from the rod in such a manner that it 
was partially submerged in water. The shape was sv 
calculated that the increase in the buoyancy just com 
pensated for the decrease in cross section when the rod 
elongated in consequence of creep. In this test the 
strain rate remained constant until fracture occurred. 
Andrade’s test suggests that Eq. (1) may suffice to 
predict the creep deformations even during the tertiary 
phase of creep if the effects of the large deformations 
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we duly taken into account.* The engineering defi 


nition of stress-namely, 0 = P/A—where P is the 
tensile load and A is the original, or nominal, cross 


sectional area, must be replaced by 


where € is the engineering strain—that is, the increment 
in length divided by the initial length. The length of 
the bar after stretching is (1 + e) times the initial 
length; hence the cross-sectional area must decrease 
by a factor (1 + e) if the volume is to remain constant. 
But such a decrease in area results in the stress given 
by Eq. (2 when P is constant. This stress is denoted 
as the true stress. 

Similarly, the definition of strain must also be re 
vised. Referring the increment in length to the initial 
length is ambiguous when a test is interrupted and 
started again. Should one use as reference the length 
at the beginning of the first or of the second phase of 
the test? A consistent definition is possible if the small 
increment in length Lde is divided by the instantaneous 
value of the length L(1 + ¢«). The increase in the nat 


ural or true strain 1s then 

de, = de/(1 + «€ (Ss 
which becomes upon integration 

€, = log (1 + e€) (4) 


This is the definition of the natural strain which is the 
natural logarithm of unity plus the engineering strain. 
When the engineering stress o and the engineering 
strain «, which are perfectly satisfactory quantities as 
long as the deformations are small, are replaced by the 
true stress o, [Eq. (2)] and the natural strain e, {Eq. 
{)] in the empirical creep law of Eq. (1), a relationship 
is obtained which should be capable of describing the 
phenomena accompanying large creep deformations. 
The relationship is expressed in the form of a differen 
tial equation. The differential equation was inte- 
grated in reference 15; from the integral the change in 
cross-sectional area can be calculated for any time. 
The results show that the cross-sectional area de 
creases to zero only after an indefinitely long time 


When 
n is greater than unity, the time necessary to reach a 


when the exponent in the creep law is unity. 


vanishingly small area is finite; it is designated as the 
critical time f,, and is given by the simple theoretical 


equation 
tery = 1/€n (5 


Here €, is the strain rate at the beginning of the creep 
test, when the engineering stress in Eq. (1) is the same 


as the true stress. 


* The following is a quotation from Andrade’s 1910 paper 
It was found that taking the rate at different times in the vis 
cous flow, and dividing by the mean length pertaining to each 
Tate, a rate per unit length was obtained which was constant at 
different stages of the flow. It was thus established that under 
constant stress the rate per unit length, once the initial effect has 
died out, is constant right up to breaking.”’ 
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STRUCTURES 


It is worth noting that the critical time docs not de 
pend on the material constant A, but only on the ex 
ponent # and the initial creep rate €&. It can be taken 
as the time necessary for fracture to occur in a tensile 
creep test. In reality, the rod fractures before its cross 
sectional area is reduced to zero; but the reduction in 
area occurs so slowly at the beginning and is acceler 
ated to such an extent at the end of the experiment 
that the theoretical difference between time elapsing 
before failure and that necessary to reach a vanishing 
area is negligibly small. Comparison with the test 
results of Dorn and Tietz gave good agreement be 
tween the theoretical value ¢,, and the experimentally 
determined time to failure. This is illustrated in Fig 
3. 
The nonlinearity of the creep law |Eq. (1)] is the 
When in 


consequence of inaccuracies of machining some section 


cause of another interesting phenomenon. 


of the rod is slightly smaller than the other sections, the 
minimum section is subjected to a slightly higher 
stress and thus experiences a materially increased rate 
of creep. Thus it contracts more rapidly than other 
sections, with the result that initial deviations from 
perfect uniformity are gravely exaggerated during the 
creep process. As no rod can be manufactured so as 
to be perfect, creep fracture always occurs with a neck 
ing down of the weakest section. 

This phenomenon was also analyzed in reference 15 
When the middle section of a rod was assumed to be 
| per cent smaller in area than the end sections and 
the variation in section along the rod was represented 
by a sine law, the progress of necking was found to be 
as shown in Fig. 4. The shape calculated closely re 
sembles those observed in experiment.® 


THE ELastic ANALOG 


The calculation of the stresses in a structure whose 
material deforms in creep is facilitated through use of 
the elastic analog. The analog was originally stated 
by Alfrey in 1944" and was extended by 
1950;'* both formulations were valid only when the 


Tsien in 


stress, the strain rate, and their derivatives were con 


nected by linear relationships. Unfortunately, in 


metals steady creep is governed by Eq. (1), which is 


highly nonlinear because m is usually greater than 3. 
It was desirable therefore to derive an elastic analog 
for nonlinear creep also. 

This was accomplished by the writer in reference 19 
for the case when the deformations of the material were 
The analog states 


completely governed by Eq. (1). 


that the stresses are the same in a nonlinear perfectly 


elastic body whose stress-strain law is 
«= (o/i)" (6 


as in the body whose deformations are completely de 
fined by the creep law of Eq. (1) provided the boundary 
conditions are suitably chosen. It is to be noted that 
the second body deforms only in consequence of creep; 
it must not have any instantaneous elastic or plastic 


deformations. The boundary conditions are easy to 
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determine. Rigid supports and prescribed loads must 
be the same for the structure that creeps and for its 
elastic analog. A prescribed displacement of a point 
on the boundary of the elastic model corresponds to a 
prescribed velocity of displacement of the same magni- 
tude in the body governed by Eq. (1). 

The analog contains no mention of strain; as a matter 
of fact, the strains have definite values, constant with 
time, in the nonlinear elastic body while they keep 
increasing with time in the body subject to creep. 

The existence of the analog was proved in reference 
19. The proof is based on the fact that in the formu- 
lation of the equilibrium and compatibility conditions 
time appears homogeneously; this means that the 
equations do not change when the time scale is changed. 

The advantage of the analog is that it permits the use 
of methods well established in the theory of elasticity 
in the solution of creep problems. Of course, the solu- 
tion of any but the simplest problems is still difficult 
because of the nonlinearity of the stress-strain rela- 
tionship given in Eq. (6). When direct integration 
of the differential equations of the problem is too diffi- 
cult, use can be made of the complementary energy 
principle. If the supports of the body are rigid, the 
principle states that the stress distribution must be 
such as to store the least amount of complementary 
energy in the system. The complementary energy 
L’ per unit volume is defined as 


ae } eda (7) 


Substitution of ¢« from Eq. (6) and integration yield 
U’ = [A/(n + 1)] (@/A)**! (8) 


Equation (6) reduces to Hooke’s law when n = 1. 
In such a case the constant \ is Young’s modulus and 


the equation can be written as 
e = (cE) (6a 
Similarly l’’ becomes the strain energy LU’: 
Ll’ = U = (07/2E) (Sa) 


The analog for the material that is subject to linear 
creep is the ordinary linear elastic material. This 
was shown earlier by Alfrey and Tsien." 

Two interesting consequences of the elastic analog 
mentioned here. First, the complementary 
energy of Eq. (8) rises much more rapidly with increas- 
Hence 


Can be 


ing stress than the strain energy of Eq. (Sa). 
the requirement that the complementary energy be a 
minimum leads to a much more radical evening out of 
the stresses in structures whose material creeps than 
the requirement that the strain energy be a minimum 
the least work principle) causes in elastic bodies fol- 
lowing Hooke’s law. In the presence of nonlinear 
creep the stresses are more uniform and the stress 
concentrations are less severe than in bodies whose 
stress-strain relationship is Hooke’s law. 

The second consequence is that an initial lack of fit 
of a redundant structure does not have any effect upon 
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the stress distribution. A prescribed displacement of 
a point on the surface of the model whose material 
follows Eq. (6) can give rise to very serious stresses; 
but its equivalent in the body whose creep is governed 
by Eq. (1) is a prescribed velocity. Hence, a constant 
rate of change of the displacement causes stresses in the 
presence of creep but not a prescribed displacement 
This is helpful in experimental work because inaccu 
racies in the test setup have little influence on the re- 
sults when the material behaves in accordance with 


Eq. (1 


THE BEHAVIOR OF REAL MATERIALS 


The statement that a lack of fit cannot cause stresses 
is startling at first hearing; it does not seem to agree 
with one’s personal experience. Of course, the state- 
ment is true only for a material whose deformations 
are caused exclusively by creep in accordance with the 
creep law of Eq. (1). Such materials do not exist in 
reality. 

At the moment of load application no creep has time 
to develop; the instantaneous deformations are gov- 
erned by the laws of elasticity and plasticity. A short 
time after the load is applied, the influence of instan- 
taneous elasticity and plasticity is still noticeable; and 
the effect of transient creep is also likely to be greater 
than that of the steady creep governed by Eq. (1). 

The theory of creep based upon the elastic analog 
or, what is the same, on Eq. (1)—1s a limit theory that is 
never completely valid in a mathematical sense but 
becomes more and more accurate as the secondary 
creep deformations progress. For engineering pur 
poses the accuracy of this limit theory may be satis- 
factory after a very short time. In reference 19 an 
example in which both the elastic and the secondary 
deformations were taken into account was 
In the example the stress distribution 


creep 
worked out. 
of the limit theory was reached within engineering 
accuracy in about one minute; at that time the creep 
deformation was about the same as the elastic defor- 
mation. 

When the stress distribution is needed shortly after 
the load is applied, more complete stress-strain laws 
must be used than Eq. (1). If both the elastic de- 
formations and the steady creep are considered in a 
a uniaxial state of stress, the law governing the de- 
formations 1s 


€=(¢6 FE) + (aA (9 


At high values of stress and temperature the elastic 
part of the instantaneous deformation may be negli- 
gibly small as compared to the instantaneous plastic 
deformation. In such a case Odqvist’s strain rate law 
of reference 20 can be written in a slightly modified 
form as 


nm 


€ = (ao/p)"(o/ nw) + (a/d)" (10 


that is, during 
When o changes 


This equation holds only for loading 
the time when o does not change sign. 
sign and thus the element is unloaded, the creep law 1s 
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Fic. 7. 


again that of Eq. (1). The equation can represent the 
instantaneous plastic deformations rather accurately 
and the primary creep deformations in an approximate 
manner. 

The three laws represented by Eqs. (1), (9), and (10) 
are compared in Fig. 5. The values of the constants 
\, w, m, and n were so chosen as to obtain the best pos- 
sible agreement with the experimental curve. If the 
period of loading had been extended to 5/4, the differ- 
ences in total strain shown in the five diagrams would 
have become negligibly small. 

Of course, it is possible to combine Eqs. (9) and (10) 
and to take into account elastic, plastic, and secondary 
creep deformations simultaneously. This was done 
in the section on creep buckling. 

It is also possible to talk of a limiting case of the 
states of creep deformation. This is approached as 
the exponent m in the creep law of Eq. (1) increases 
In the analogous elastic body whose 
(6), » changes in the same 


beyond all limits. 
deformations follow Eq. 
manner. 

Stress-strain curves of the analogous elastic material 
are shown in Fig. 6. It can be seen from the figure 
that the behavior of the material approaches that of a 
rigid-plastic, or perfectly plastic, material when » in- 
creases beyond all limits. Hence in the limiting case 
of all the steady creep states, the stresses are distrib- 
uted in the body in the same manner as they are in a 
perfectly plastic body. This limiting case can there- 
fore be studied by means of the well-developed methods 
of perfect plasticity and limit analysis.”! 


STRESSES IN A PIN-JOINTED FRAMEWORK 


Fig. 7 is a schematic representation of a pin-jointed 
framework on which a load P is acting. The stresses 
in this framework will now be calculated because 
frameworks are the simplest structures to analyze. 
It will be assumed that the material of the framework 
follows the creep law of Eq. (1), and use will be made 
of the elastic analog. 

If F stands for the tensile force in the slanting bars 
and X for that in the vertical bar, equilibrium of forces 
requires that 


2Fcosa+X =P (11) 
that is, 


F = (P — X)/(2 cos a) (12) 
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The magnitude of Y can be determined from thy 
complementary energy principle. 
that the length of the slanting bars is (/,/cos a), th 
total complementary energy UL’ stored in the thre 


bars is from Eq. (8) 


If it is remembered 


U’ = [hy/(n + 1)] jAx(X/AxdA)"*! + 
(2A ~/cos a) [((P — X)/2A rd cos al"t!} (13 
where Ay and A, are the cross-sectional areas of the 
vertical and the slanting bars, respectively. 
In consequence of the complementary energy prin 
ciple the derivative of U’’ with respect to Y must vanish 
ol’ /oxX = 0 14 
Differentiation and solution for X yield the expressior 
: 2 (1 + 2(A, Ax) (cos a)” te (15 
When the cross-sectional areas of the bars are equal 
this expression reduces to 
x = P/(1 + 2Aecosa)”""") (15a 
Fig. S shows the variation of the ratio x of the force ir 


the vertical bar to the applied load 


x = (X/P) 16 





for a framework in which a = 60° when the exponent 


n in the creep law changes. The linear creep law, } 
n = 1, corresponds to the case of the perfectly elastic | 
material for which Hooke’s law is valid. In a frame- 
work built of such a material, the vertical bar carries 
SO per cent of the load. 

As n increases, the share of the vertical bar is re- 
duced. This is in agreement with the observation 
made in the preceding section. When ” becomes very 
large, x approaches 0.5. Then half the load is carried 
by the vertical bar and the other half by the two slant- 
ing bars. But the force F in the slanting bars is the 
same as that in the vertical bar, because the vertical 
component of Fis 0.5/. This force distribution repre- 
sents the utmost in leveling out the loads in the bars 
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of the framework. It corresponds to the stress dis- 
tribution assumed in limit analysis. 

Next the effect of the elasticity of the material, 
which was disregarded in the analysis jiist completed, 
will be considered. Eq. (1) will be replaced by Ea. (9). 
The condition of equilibrium [Eq. (11)] remains un- 
changed. The condition of the compatibility of the 
deformations can be written directly; it replaces the 
principle of the minimum of the complementary energy, 
which is not valid for the law of deformations expressed 
in Eq. (9). 

The elongation of the vertical bar can be designated 
is Ay; it gives rise to a vertical downward displacement 
Ay of the point of application of the load P [see part 
b) of Fig. 7] from the original position A to the final 
position C. The elongation A, of the slanting bar must 
be such magnitude that the continuity of the frame- 
work is maintained. During the deformations the 
slanting bars rotate about their upper pin joints; as the 
angle of rotation is small, the path followed by the 
point of application of P can be taken as perpendicular 
to the original direction of the slanting bar. Without 
the rotation the lower end point of one of the slanting 
bars would have arrived at point B; after the rotation 
it must be at C. The condition then is that 


Ar/Ax = cosa (ig 
Sut 
Ay = [és h dt (1s 
J0 
ind 
Ar = ec cos ajdt (19 
J0 


As. Eq. (17) must hold at all times, one obtains 
€r/ €x = COS’ a (20 


Equations (9), (12), and (20) can be solved simul- 
taneously. One obtains for the time ¢ necessary to 
reach a state in which the vertical bar carries x times 


the applied load P 














Fic. 9. Framework tested at 600° F. 
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Fic. 10. Load history in vertical bar of pin-jointed framework 
P = 400 Ibs 
my dx 
t=A | (21 
Jz, (1 — x)" — Bx’ 
where 
A = ${2"(cos a)"*? [1 + (1, 2 cos® a)]}/[(E/a)(o/d)"] 
(2la 
B = 2"(cos a)"*t? (21b 
» «i Z 
, {.1ic 
c= X/P) 


and x,; is the value of x corresponding to the linear 
elastic distribution of the load and xy,;, that corre- 
sponding to the final distribution for which ¢ is sought. 

The integration can be carried out graphically 
without difficulty and x can be plotted against time. 
The initial value of x when ¢ = 0 corresponds to the 
linear elastic distribution; at the moment when the 
load is applied Hooke’s law is valid because it takes 
time for the creep deformations to develop. When 
the denominator in the right-hand member of Eq. (21 
vanishes, ¢ increases beyond all bounds. This happens 
at 


; 1/7 2» 
“limit 1/1+B } \<< 


Consequently, x approaches asymptotically the 
value Xjjmit aS time increases beyond all limits. It is 
easy to show that this value is the same as the one 
corresponding to Eqs. (15a) and (16). Thus the cal 
culations carried out by means of the elastic analog 
on the basis of Eq. (1) represent the asymptotic solu- 
tion of the more complex problem in which the elastic 
deformations, as well as the secondary creep deforma 
tions, are considered. 

Several frameworks have been constructed and 
tested at the Polytechnic Institute of Brooklyn to check 
the validity of the theory just presented (see reference 
22). Some of these frameworks had rigid joints while 
others were provided with pinned attachments. The 
material of the frameworks was 24S-T4 sheet of 0.064 
in. thickness. The width of each bar was '/» in. and 
h was 10in. 

The test setup is shown in Fig. 9. 
was taken after the vertical bar necked and broke. 


The photograph 
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Fic. 11. Creep bending test. 


Before the actual test was begun, the pin-jointed frame- 
work was placed in an oven and heated by a number of 
electric resistance heating elements. The upper ends 
of the slanting bars were attached to adjustable sup- 
ports, and the vertical bar was suspended from a 
dynamometer. The temperature of both oven and 
specimen was raised to 600°F. in a heating period of 
about an hour. When the thermocouples indicated 
a uniform and steady temperature, a weight of 400 
Ibs. was suspended from the bottom joint. Readings 
were made of the temperatures during the test process 
to insure a constant temperature of 600°F., and the 
load indicated by the dynamometer was recorded. 
Results obtained from one such test are plotted in Fig. 
10. 

In order to compare the predictions of the theory with 
experiment, the integral in Eq. (21) had to be evalu- 
ated for the framework tested. The material constants 
needed in the calculation were obtained from bending 
tests. Fig. 11 is a photograph of the test setup. 

An aluminum alloy beam of 12-in. span was supported 
on knife edges and was loaded by two symmetric ver- 
tical loads in such a manner as to produce constant 
bending moments between the supports. Changes in 
the curvature of the beam were measured by means of 
heat-resistant electric wire gages (see reference 23) 
bonded to the faces of the beam. The beam was 
placed in an oven whose temperature was maintained 
uniformly at 600°F. The tests indicated that the 
material constants of the beam were 


9 


n=3 t (23) 
\ = 87,500 Ibs. in.-2 hr.'/3 f _ 


The elastic modulus of the material was determined 


from compression tests made with columns and was 
found to be 


E = 7.4 X 108 lbs. per sq.in. (24) 


With these values theory yielded the dotted curve of 
Fig. 10. 

Comparison of the theoretical and experimental 
curves showed satisfactory agreement only during the 
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first third of the test period. The initial SO per cent 
load in the vertical bar dropped relatively rapidly 
during the first hour of testing. 
in the test became smaller afterward, but not so smalj 


The asymptotic value of the 


The rate of decrease 


as predicted by theory. 
theory, x = 0.612, was reached in the experiment after 
2 hrs., 15 min., and x decreased further without appar- 
ently approaching any asymptote. After about § 
hrs. of testing, the load started to decrease at an ac- 
celerated rate and the vertical bar fractured. 

There is not enough information available at present 
to explain the deviations between theory and experi- 
ment. It must be remembered, however, that the vyal- 
ues of the material constants were not too reliable be- 
cause they were obtained from beam and column tests 
carried out with drawn sections while the bars of the 
framework were made of sheet material. Moreover, 
the history of the heat treatment of the material dif- 
fered; also 24S-T3 material is somewhat unstable and it 
may have changed during the 8 hrs. of testing. Un- 
doubtedly, the substantial drop in the load carried by 
the vertical bar during the final stages of the experi- 
ment was caused by the large elongations accompany- 
ing necking. This necking, and the subsequent frac- 
ture, were predicted by theory for a much later period 
of time. It is possible that necking was precipitated by 
changes in the structure of the material as well as by 
inaccuracies in the manufacture of the specimen. 


CREEP BUCKLING 


A discussion of creep buckling must begin with the 
observation that every column is imperfect. Its center- 
line is straight only within the accuracy of the machin- 
ing operations, its material is inhomogeneous, and its 
alignment in the testing machine or its attachment to 
other elements of a structure can never provide per- 
fectly central load application. If it were possible to 
establish perfect conditions, the column would not 
buckle in spite of the creep deformations. Its cross 
sections would increase and its length would decrease 
with time under the action of the constant axial com- 
pressive load, but the perfect symmetry of the system 
would always be preserved. The column would finally 
resemble a pancake in shape, but only after an indef- 
nitely iong time would its length be reduced to zero. 
that is, imperfect— column 
The distance y be- 


The behavior of a real 
is entirely different (see Fig. 12). 
tween the line of action of the constant compressive 
force P and the centroid of a cross section of the column 
serves as a lever arm which multiplied by the force is 
the bending moment acting on the cross section. This 
moment causes creep when the temperature of the 
column is high enough; the creep deformations, in turn, 
increase the initial deviations from straightness. When 
the exponent m in the creep law [Eq. (1)] is large, the 
column behaves in a typically nonlinear manner. Its 
curvature increases very slowly at the beginning of the 
loading. The lateral deflections become noticeable 
after a while, and soon after the column snaps through 


and collapses. 
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Perhaps the most interesting feature of creep buckling 
is the absence of a critical load. Creep deformations in 
creasing the initial deflections take place continuously 
ynder any compressive load however small. Thus the 
column buckles eventually even if the loading acting 
on it is much smaller than the Euler buckling load cal- 
culated from the elastic properties of the material at 
the particular temperature. The important quantity 
is the length of time which elapses between the appli- 
cation of the load and collapse. The designer must 
make sure that this time, designated as the critical 
time, is greater than the lifetime of the structure of 
which the column forms part. 

Creep buckling was first analyzed by Freudenthal 
in 1946:24 at about the same time it was also investi- 
gated experimentally by Ross.” In the United States 
the first known experiments in this field were carried 
out by Marin in 1947," and investigations at the 
Batelle Memorial Institute (1949)? and the National 
Advisory Committee for Aeronautics followed suit 
1952). The author and his group at the Polytech- 
nic Institute of Brooklyn took up the problem of creep 
buckling in 1949;°° the investigations were continued 
under a project sponsored by the National Advisory 
Committee for Aeronautics, and various aspects of the 
problem were clarified by Kempner, Patel, Pohle, and 
others in a series of reports.*” 

In his Wilbur Wright Memorial Lecture, the author?! 
presented a theory of creep buckling based on the 
issumption that all the deformations could be attri- 
buted to steady creep governed by the law of Eq. (1). 
Kempner’s work® has shown, however, that the elas- 
ticity of the material is capable of modifying the de- 
formations of the column sufficiently to result in con- 
siderable differences in the predicted critical time 
when the load is near the Euler load. Moreover, 
Odqvist** has pointed out that the author’s equations 
in reference 31 could easily be integrated even if the 
instantaneous plastic deformations were included in 
the analysis and that the transient part of the creep 
deformations could also be considered in an approxi- 
mate manner. On the basis of this information the 
author developed a new theory in which the various 
effects mentioned were taken into account simul- 
taneously.** The most important results of this theory 
will be presented and discussed here. 

First the stress-strain rate relationship was assumed 


in the form 
é = (¢/E) = (a/p)™ (c/n) & (c/d)* (25) 


This relationship can be recognized as a combination 
of Eqs. (9) and (10). The first term in Eq. (25) is 
clearly the contribution of the elasticity of the material 
at the test temperature. The remaining constants in 
the expression can be determined without difficulty 
lor any material for which reliable creep diagrams of 
the type of Fig. 1 are available for a number of constant 
loads or stresses. 

If in Fig. | the straight line portion of the creep curve 
is extended until it intersects the elongation axis, the 
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intercept on the elongation (or strain) axis represents 
the instantaneous deformations following the load 
application. Integration of the first two terms in the 
right-hand member of Eq. (25) gives a relationship 
which is valid at ¢ = 0. 

e = (¢0/E) = [1/(m + 1)] (o/p)* +! (26) 
The sign in the equation must be suitably chosen for 
each condition of loading and for each value of m. 
For instance, the negative sign is needed when a com 
pressive (negative) stress is applied and m = 1; in such 
a manner the terms representing the elastic and the 
plastic parts of the instantaneous deformations both 
make negative—that is, compressive—contributions 
to the strain e. 

The elastic part of the intercept on the elongation 
axis can be computed as soon as £ is known for the par- 
ticular temperature. The remainder must then be 
due to the second term in the right-hand member of 
Eq. (26). One can subtract (¢//) from the intercept 
and prepare a double logarithmic plot of the differ- 
ence in strain against the stress. If the points rep- 
resenting all the creep tests can be connected by a 
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Fic. 13. Static loading of slender column 


straight line, the slope of the line is m + 1. This 
graphic procedure is the same as the one used to de- 
termine the value of m from minimum creep rate data. 

In reference 32 Odqvist assumed n = 3 and m = 1. 
On the basis of tests now in progress at the Polytechnic 
Institute of Brooklyn, these values seem to represent 
approximately the deformations of 24S-T3 aluminum 
alloy bars at 600°F. In addition the yet incomplete 
data indicate that the choice 


A = $7,500 
800,000 (27) 
E=74*X 104 


= 
II 


is a reasonable one when the constants are measured in 
the units inch, pound, and hour. 

A detailed derivation of the equations of the static 
loading and of the subsequent creep process can be 
Unfortunately the irreversible 


found in reference 33. 
plastic deformations necessitate many changes in the 
equations and thus lead to an inconveniently large 
number of final expressions, each of which is valid only 
during a particular phase of the loading or creep 
process. Instead of reproducing these equations here, 
the writer wishes to present a numerical example that 
was completely worked out for a column of a slender- 
ness ratio of 111. For the sake of simplicity the col- 
umn was assumed to have two concentrated flanges 
each of a cross-sectional area A/2. The flanges were 
kept a distance h apart by a perfectly shear-resistant 
web of no area. This highly idealized column was 
shown to represent actual conditions with a satisfac- 
tory degree of accuracy. 

The initial shape of deviations was assumed to be of 
the half sine wave type with an amplitude of 0.06 p, 
where p is the radius of gyration of the cross section of 
the column. The ratio of the amplitude to p was des- 
ignated by the symbol dg. It was shown in reference 
31 that the half sine wave shape was a statisfactory 
approximation to the actual conditions in the analysis 
of creep buckling. 


As far as the instantaneous response to the load a 
plication is concerned, one reservation must be mag 
Even though the load is applied rapidly, it must no; 
be allowed to give rise to dynamic effects. In an e 
periment, the supports of the load must be remoy, 
carefully to prevent impact and vibrations. From th, 
standpoint of the theory, it is required that the tin, 
interval between the no-load and the full-load cond; 
tions must be so short that the effect of the thir 
term in the right-hand member of Eq. (25) should } 
negligible. 

Fig. 13 presents the results of the analysis of the 0 
called static loading process which satisfies the cond; 
tions just stipulated. During the first phase of th 


loading the compressive stress increases in both flanges 


of the column. All the deformations are governed 
therefore by Eq. (26) with the negative sign before th 
second term because m = 1. When a nondimensiona] 
deflection @ = 0.241 is reached, the compressive strain 


in the flange on the convex side of the column begins 
to decrease. The point corresponding to the strain 
reversal is indicated in the diagram. 

When the compression decreases, only the elasti 
part of the deformations can be regained. Conse 
quently, the deformations in the flange on the conve) 


side are governed by the equation 

e = (0/E) — (1/2) (¢,/p)? (28 
where 

o, = —4,242 lbs. per sq.in. (28a 


is the average stress prevailing at the moment of strair 
reversal, while the stress-strain relationship for the 
flange on the concave side remains Eq. (26) with m = 
1. These two relationships define the deflections of 
the column if the static loading is continued. A 
second change in the equations is required when the 
compressive stress in the flange on the convex side 1s 
reduced to zero. From then on the stress-strain law is 


e = (a/E) — (1/2) (o,/p)? + (1/2) (o/p {29 


In the numerical example tension begins to develop at 
the middle of the convex flange when the nondimen 


sional deflection is 
a& = 1.0 (29a 
and the average stress 1s 
o = (P/A) = —5,285 lbs. per sq.in. (29b 


Unfortunately, the actual conditions are even mort 
complex because the average stress must be further 
increased before tension can spread from the middle 
of the convex flange toward the end points of the 
column. If this is disregarded and if Eq. (29) is as 
sumed to be valid for the entire length of the convex 
flange as soon as the average compressive stress exceeds 
the absolute value of 5,285 Ibs. per sq.in., the remainder 
of the curve in Fig. 13 is obtained. 

This curve has a maximum of —¢ = 5,381 Ibs. pet 
sq.in. at d@ = 1.922. The existence of a maximum 
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indicates instability. If the load is increased beyond 
the maximum, no static equilibrium is possible for the 
column. The excess load beyond the resisting capacity 
of the column causes lateral accelerations, and the 
column buckles suddenly, usually with an audible thud. 
Hence for static loading the critical stress in the column 
of the numerical example is 


Scr = —9,381 lbs. per sq.in. (30) 


0 


The creep behavior of the column was next investi- 
rated. It was assumed that the static loading was 
halted at an average stress 

o¢ = P’A = —4,500 Ibs. per sq.in. (31) 
[he corresponding nondimensional static deflection 
This value is beyond that 


was found to be dy = 0.298. 
at which strain reversal occurs but below the value at 
which tension begins in the flange on the convex side. 
The load causing the average stress given in Eq. (31) 
is 83.5 per cent of the static buckling load. It was 
assumed to be maintained unchanged, and the column 
was allowed to develop increasing deflections in con- 
sequence of creep. 

The creep law of Eq. (25) was taken as the basic re- 
lationship governing the lateral deflections of the col- 
umn. In each phase of the creep process the signs had 
to be chosen in a suitable manner, and the term repre- 
senting the instantaneous plastic deformations had to 
be omitted whenever a flange was being unloaded rather 
than loaded. A detailed analysis of the various phases 
of the creep process will not be given here, but the re- 
sults of the calculation are presented in Fig. 14. 

It can be seen from the diagram that the column 
buckles in about 43 min. after the load application. 
In the analysis the critical condition is characterized 
by an infinite velocity of the lateral deflections. Be- 
cause of the many changes in the relationships govern- 
ing the deformations, it is not possible to give a single 
that is, the 
time at which the velocity becomes indefinitely large 
can be computed. It is likely, however, that the formula 


expression from which the critical time 


Aer = 2.36 A/(o/p)? (32) 


gives a reasonable, although somewhat unconservative, 
estimate of the critical nondimensional deflection when 


m=landn= 3. In this equation 


A = ex + (¢/E) — (1/2) (o/u)? (3333) 

and 
€x = m?/(L/p)? (34) 
is the Euler strain and ¢ = P/A is the average stress, 


considered negative when compressive. The critical 
time elapsing between load application and collapse can 


then be computed from the equation 


a —=( 7/3) JA, er? 4 + Apo’ 
(o/d)* (2 we ay” 4 + a,’ 


s) nd ) 
2B tan 1 | 2Ger a) i (35) 
[ $ + d,,dy if 
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with 

B = (4/32) (o/p)? (35a 
where log indicates the natural logarithm. In this 
equation again d is the nondimensional deflection 
amplitude of the column immediately after the load 
application before the steady creep deformations be 
come noticeable. As a compressive stress is considered 
negative, Eq. (35) yields a positive critical time. 

In Fig. 15 the effect of the magnitude of the applied 
ioad upon the critical time is shown for the particular 
column investigated. The curve labeled /,, , represents 
the upper limit on the critical time computed from Eqs 
(32) to (35). The lower limit is designated by f 
it was calculated from the theory of reference 35 
The third curve, labeled 4),,;¢, indicates corresponding 
values of the time and average applied stress at which 
tension begins in the flange on the convex side of the 
column. 

Probably the most interesting feature of the creep 
buckling process revealed by Fig. 15 is the overwhelming 
influence of the average applied stress upon the life 
time of the column. A slight increase in the cross 
sectional area of the column can therefore contribute 
materially to the safety of the structure. The close 
agreement between the upper limit /,, , on the critical 
time computed from the simple equations given earlier 
and the lower limit /,, ; should not be taken for granted. 
The deviation between the two curves becomes greater 
when the slenderness ratio (ZL ‘p) of the column and the 
value of the material constant yu are decreased. 
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A Critique of Shock-Expansion Theoryw’ 


J. J. MAHONY?T 
University of Manchester, England 


SUMMARY 


The ptions underlying the application of shock-expansion 


rv to the calculation of the pressure distribution on a thin 


two-dimensional airfoil in a supersonic stream ar 


rp-nosed 


is suggested that they might be expected to lead 


mined, ana tt 
ippreciable errors when the Mach Number is large rhis 
is then examined by the use of the “‘hypersonic analogy,” 
{it is shown that for circular-are airfoils the pressure distribu 
is given to good accuracy by shock-expansion theory even 
1 the quantities neglected are no longer small An asymp 
otic form for the decay of the leading-edge shock is developed 
the case of shocks that are too strong initially for the Fried 
hs theory to apply 


(1 INTRODUCTION 


| PRESSURE DISTRIBUTION on a two-dimensional 
thin airfoil with a pointed nose in a uniform super- 
sonic stream is normally calculated by the use 
In this theory, it is assumed 


of 
shock-expansion theory. 
that the shock waves attached to the leading edge are 

f uniform strength in that portion of the flow field 
which affects the pressure distribution on the airfoil. 
Furthermore, it is assumed that all boundary-layer 
effects may be ignored. Then the flow pattern in the 
neighborhood of the airfoil becomes that of a homen- 
tropic simple wave, and, for a given section, the pres- 
sure distribution is immediately deducible. 

Both of these assumptions are in general only ap- 
proximately true. The boundary-layer effects will, 
i course, be small apart from near the trailing edge 
where the boundary layer may be separated due to 
The 
eflect of this on the pressure distribution cannot be 
calculated analytically, although a fair amount of ex- 
perimental information (for example, references | and 


the upstream influence of the rear shock wave. 


2) exists from which the rise in pressure near the trailing 
edge can be estimated in many cases. The assumption 
regarding the leading-edge shocks will, in fact, be true 
lor a double wedge section for sufficiently large Mach 
Numbers, but for all sections in which the surface is 
curved near the leading edge, the shock will decay 
within the region of interest due to the influence of the 
expansion wave generated by the airfoil. Thus, shock- 
expansion theory might be expected to yield accurate 
results only if the decay of the shock and the conse- 
juent rotationality of the flow fleld could be neglected. 
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Thus, one is led to inquire whether, under any condi 
tions of thickness ratio and free-stream Mach Number 
liable to be encountered in practice, the errors may 
become appreciable. 

Such experimental investigations*~* as have been 
carried out indicate clearly that, for conventional thick 
ness ratios of high-speed sections and for the Mach 
Numbers of these tests, the theory gives satisfactory 
results for lift and wave drag. However, the investi 
gations that have been made do not exhaust all possible 
combinations of Mach Number 


particularly at extremely high Mach Numbers 


section, 

Now, 
order of magnitude arguments suggest that for a given 
thickness ratio, the with Mach 
Number; for the changes across a shock wave are not so 


and airfoil 


errors will increase 


much functions of the stream deflection, at least as 
far as the usual Busemann expansions are concerned, 
as functions of the product of this with the Mach Num 
ber. This product could possibly be as large as unity 
for a thin sharp-nosed airfoil, and yet the Mach Number 
may be within the range where the gas may be treated 
as a continuous medium. Thus, in the case of high 
Mach Numbers, it is not possible to use the fact that 
the entropy change is of third order in the stream de 
flection to deduce that the effects of the vorticity be 
hind the shock are small. 

Thus, it seems that a further investigation of the 
case of high Mach Numbers is necessary. Further 
more, it would seem reasonable that if the rotational 
nature of the flow can be neglected in this more crucial 
case, then it can be neglected also at lower Mach Num 
Thus, as a preliminary investigation, attention 
The mathe 


bers. 
will be confined to the hypersonic case. 
matical difficulties arising in the treatment of a rota 
tional supersonic flow may be reduced considerably if 
an appeal is made to the hypersonic analogy due to 
Hayes? Goldsworthy. In_ this 


theory an analogy is drawn between the airfoil prob- 


and developed by 
lem and the one-dimensional motion of a perfect gas 
due to a certain calculable motion of a piston that 
bounds the gas on one side. The physical basis for this 
analogy is that, because the Mach angles are small, 
velocity changes and gradients in the direction of flow 
are smaller than those in the perpendicular direction. 
This implies that planes of fluid initially normal to the 
stream move essentially in their own plane under the 
laws of one-dimensional unsteady motion as they move 
downstream. For the pointed-nose airfoil under con- 
sideration, the analogous unsteady problem is that of 
a piston set impulsively into motion with a velocity 
whose ratio to the undisturbed sound speed in the 
gas is thereafter equal to the product of the Mach 
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Number and the stream deflection in the original prob- 
lem. The main part of the paper is devoted to a 
study of this unsteady problem. 

The changes across a shock wave as functions of the 
shock strength are displayed in Table | for values of 
y, the ratio of specific heats of the gas, equal to 1, 7/5, 
and 5/3. It is readily apparent that the increase in 
entropy and the change in the Riemann invariant, 
whose variation is neglected in shock-expansion theory, 
remain small even for the comparatively strong shocks 
involved. They are, in fact, much smaller than the 
third-order Busemann expansion would suggest. Thus, 
the flow will not differ greatly from an irrotational 
simple wave. However, as is shown in the next sec- 
tion, the pressure distribution is rather sensitive to 
small variations of this Riemann invariant. Thus, 
complete neglect of these variations would be unjusti- 
fied, but a solution in which their squares were neglected 
would appear satisfactory. Such a solution can be 
obtained if the shock wave is calculated from a simple 
wave theory and the resultant entropy is distributed 
along streamlines calculated in a like manner. Then 
a second approximation to the flow pattern can be ob- 
tained with this calculated entropy distribution and 
initial flow corresponding to the calculated shock. In 
theory, this process could be repeated to obtain higher 
order approximations, but this will be unnecessary for 
the present purpose. 

The pressure on the piston at a given instant is known 
if the value of the isentropic Riemann invariant is 
known. This has a small variation along the charac- 
teristic, which reaches the point from the shock wave 
owing to the rotationality of the flow field. The value 
at the shock wave is also different from that assumed 
in shock-expansion theory. It will be found that the 
numerical values of the departure from shock-expansion 
theory at the shock wave and the change along the char- 
acteristic cancel out almost exactly in the examples of 
circular-are airfoils. It thus appears that the accuracy 
of shock-expansion theory in the extreme cases, when 
the product of the Mach Number and the maximum 
stream deflection is as large as unity, is to some extent 
a matter of chance. In a final section the influence of 
finite Mach Number is discussed briefly. 

The general conclusion of the paper is that shock- 
expansion theory may be used with confidence to pre- 
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dict the inviscid flow about an airfoil in a hypersonj 
stream even under conditions where the shock-way, 
decay, in the region that can influence the surfag, 
pressure distribution, is considerable. 


(2) MATHEMATICAL FORMULATION 


Consider the one-dimensional motion of a_ perfect 
gas initially at rest, when it has a uniform sound speei 
¢;, in a semi-infinite tube bounded at one end by a pis 
ton which at zero time begins to move with a spee 


6c. The piston is then subject to a non-zero retarda 


tion (not necessarily constant) until the piston returns 


to its initial position where it is brought to rest.7  ]j 
a is a typical value of the retardation, uc), the velocity 
of the gas, cc,, the local speed of sound, ®, the ratio oj 
the specific entropy—-with the initial state taken as 
zero level— to the specific heat at constant volume, pp 
the pressure, and pp;, the density,{ all measured at 
time ct/a, at a point whose distance along the tuh 
from the initial position of the piston is c,*x ‘a, then th 
characteristic form of the equations of continuous 


motion are: On any line 


dx/dt =u+e 
da = ,C [2¥(¥ = Ll) ]i de® f Bis 
On any line 
dx/dti = u—c i 
dB = —}c/[2y(y — 1)]{ de} ay 
and on any line 
dx/dt = ul as 
dp = OF — 
where 
a = (u/2) + c/(y — 1) (2 
and 
B = (u/2) — c/(y — 1) 2.5 


For a perfect gas the variables of state are related by 


t This is analogous to the zero incidence case in the airfoil prob- 
lem, but the restriction involves no real loss of generality as the 
extension to the case of positive incidence is straightforward 

{t Where the subscript 1 denotes conditions in the initial stat 
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2 y~i ‘ 
( exp [—@ = l (4 


The boundary conditions, subject to which the above 


ations have to be solved are: (1) on the piston path 


“= dé dt (2.8 


\ E(t), 


ind (2) on the shock wave, which is the line 


dx /dt X(6) through x = 0, ¢ = 0 2.9 
ul 6, ¢ V(6) 6 = (6 (2.10 
where 
X = sign 6, [(y + 1)/4]6 + 
Vist [(y +1) 47°84 (2.11 
} V1I+() 1) [¥6 — (1 2)62] (2.12 
ind 
® = In} V2{(X — 6) /X]~"! (2.13) 
Since the entropy is a known constant ®(4)) = 


on the piston, which follows a particle-path, until the 
piston encounters another shock, the pressure on the 
piston during its motion can be calculated from Eqs. 
9.5) and (2.7) provided 8 can be determined. In 
shock-expansion theory, 8 is assumed constant with the 
value By = B(do) 
easily be seen that the fractional error in the pressure, 


From Eqs. (2.5) and (2.7) it can 


due toa small error A@ in 8, is given by 
Ap/p = —(278/c) (48/8) 


and hence that an error of | per cent in 6 produces an 
error of the order of 7 per cent in the pressure near the 
beginning of the motion and even larger percentage 
errors near the end when y = 1.4. Thus, it is im- 
portant that the value of 8 used in calculating the pres- 
sure be quite accurate. 

The actual value of 8 on the piston can be obtained 
by integrating Eq. (2.2) along the appropriate char- 
acteristic from the shock to the piston. Hence, in the 
problem under consideration, shock-expansion theory 
has errors arising from two sources. Firstly, it applies 
the wrong initial value of 8 on the shock in that it does 
not allow for the shock decay produced by the expansion 
wave emanating from the piston. Secondly, it neglects 
the variations of entropy in the integration of Eq. 
2.2). Since the effect of the decay of the shock is 
to increase the initial value of 8 and the effect of the 
entropy variations is to decrease 6, these two sources 
of error in the determination of the value of 6 on the 
piston will tend to annul each other. Thus, one could 
expect shock-expansion theory to be more accurate 
than either of the two inherent errors would suggest 
separately. To determine the extent of this cancel- 
lation of errors, it is necessary to integrate Eq. (2.2) 
along the various characteristics. This 
the determination of the shock wave and the stream- 
line pattern, for both of which a knowledge of the dis- 


necessitates 


tribution of 8 is required. It is here that use can be 
made of the smallness of the entropy variations, for, 
if the entropy distribution 1s calculated by a homentro 
pic simple wave theory, the error involved in using 
this distribution in integrating Eq. (2.2) will be of the 
order of the square of the entropy variations and hence 
negligibly small. 

The form of the characteristic equations |Eqs. (2.1), 
2.2), and (2.3)], while convenient for numerical pur 
poses, is not particularly suited for analytical work. 
Characteristic parameters are more suited to this pur 
pose, but, because of the nonintegrability of Eqs. 
(2.1) and (2.2), no suitable parameters are available. 
In shock-expansion theory a plays such a role and hence 
can be expected to play approximately the part of a 
characteristic parameter. Since the flow is not so 
greatly dissimilar to a simple wave, 8 cannot be ex 
pected to play such a useful role, and hence another 
parameter must be found to distinguish the members 
of the other family of characteristics. For the moment 
this will be denoted by s and will be limited only by the 
condition that it is constant along, and distinguishes, 
any particular line dxdt = u — c. The definition will 
be completed later in such a way as to simplify the 
mathematics as much as possible. Then, apart from 


the possibility of the existence of limit lines and edges 


of regression,? each point between the ‘leading- and 
trailing-edge’’ shocks will be specified by the local 


value of a and s. Thus, all the flow variables and 
x and ¢ may be regarded as functions of these two vari 
ables. 

It is, however, more convenient to replace x» and f 
by a set of three dependent variables, /, G, and //, de 


fined by the relations 


dt = Fds + Hda 2.14 
dx = (F + G) (u + ¢ lds + TT(u ( da 2 ta 
where use has been made of the definition of s. It 1s 


worthy of note that the specification of F or s along any 
line on which s is not constant completes the definition 
Here, the definition of s will be completed 
Since 


of the other. 
by prescribing F along the leading shock wave. 
both dx and dt are perfect differentials it follows that 

(2.16 


oll Os = oF Oa 


and 
3) : ‘ y+ 1 a= 4 
(F +G = a + r rv) = 
a 2 Z 
3) Sa a 
i ( a+ 
Os ») ») 


where use has been made of the definitions of a and 8. 
The 


(2.1), (2.2), and (2.3) ] become 


characteristic compatability conditions [Eqs. 


08 0a = —jc/[2¥(y — 1)]} (O® 0a 2.18 


0? 0a + jLlc [Fo + (u + « G]{ (O0® Os) = 0 (2.19) 


+ These do not arise in the problem under consideration, but 
in any problem where they did they could be treated in an an 
alogous fashion to their counterparts in steady plane flow.’ 
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and 


G = 2c?H(O’/Os) + 


}(u +c) [2y(y — 1) — c(0%/da)]} (2.20) 


when written in terms of the new variables. The set 
of Eqs. (2.16) to (2.20) form a sufficient set for the de- 
termination of the five functions F, G, H/, 8, and ® 
subject to the boundary conditions on the shock wave 
and piston path. 

To obtain the boundary conditions in an appropriate 
form, one observes that the piston path dx = udt 
through x = 0, ¢ = 0 maps onto a line 


ds/da = cH/|cF + (u + c)G] through s = 0 


a=a (2.21) 


provided the origin of s is suitably chosen. Moreover, 


on this line the time is related to the velocity by 

t= T(u) = T(a + B) (2.22) 
which is the inverted form of Eq. (2.8). Hence, from 
the definition of F and // it follows that 


H = T'(a + B) [1 + (03/0a)] }F + [(u/c) + 1]G} + 
»2F + [(uc) + 1]G — T’(0B/ds)} (2.23) 
on the line that is determined by the differential equa- 
tion 
ds/da = T’(a + B) [1 + (08/0a)]+ 
}2F + [(u c) + 1]G — T’(op Os)! 


(2.24) 
and the initial condition s = 0 when a = aq. 
Similarly the differential equation governing the 
shock wave is obtained to be 
(da/d6) (dé/ds) = da/ds = 
(F/H) ((6 + ¥ — X)/(X¥ + Y — 6)] + 
(G/H) [(6 + VY) /(X + Y — 4)] 
If G, Ob Os, and Of Oa are eliminated from Eqs. 
(2.19), (2.20), (2.25), and the equation 


(2.25) 


(d@ d6) (db/ds) = d® Os = 
(OP Os) + (O®/ 0a) (da dé) (dé ds) 


for the rate of change of entropy along the shock, a 
quadratic equation for dé/ds is obtained. This equa- 


tion has the two solutions 


dé/ds = —(F/I) (dé/da) (2.26) 
and 
dé e Fo+Y-X E V*(d@ /d6) | 
ds HX+YV—8/ ld’ 2y(y — 1) (X — 8) 
= (F/I)y(6) 
(2.27) 


The first solution is rejected as physically unreal for it 
predicts that a shock wave followed by an expansion 
(F/H negative) would increase in strength due to the 
interaction, and this is contrary to experience. Fur- 
thermore, Eq. (2.26) predicts a nonvanishing value of 
dé/ds for 6 = 0, and this, in conjunction with previous 
remarks, would imply that an expansion wave would 
be shock producing. Thus, although the reason for 
the existence of this solution is obscure, the author 
feels justified in ignoring it. 
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The other solution [Eq. (2.27) ] possesses no such o} 
vious inherent physical difficulties and is, therefo,, 
the desired solution. 
in the pseudo-characteristic plane is determined } 


Thus, the equation to the shoc 


the integration of this differential equation subject 
the initial condition 6 = 6) when s = 0. Once th 
shock wave has been determined, the distribution , 
® and # throughout the flow field can be calculate 
from the streamline and characteristic pattern ay 
the characteristic compatibility conditions. For j 
stance, 8 can be determined by rewriting Eq. (2.18) ; 
the form 


0B/0a = [(8 = Ge) ty | (OP Oa) 


which is a first-order ordinary linear differential equatio; 
and can be integrated in the normal fashion. Thus 
the expression 


l 
B(a, S) = Bx(s) exp J [P(a, 5s) — Dy(s) it + 
l4y f 


{" ee $1 
Jans” EN”? V4 


where the subscript * denotes values on the shock way. 


P(a, s) — P(E, s) l dé (2.98 
f 
¥ 


is obtained for . 
If F, G, and H were known, the shock shape and th 
distribution of ® and 8 would be readily determinabk 


but since the determination of the first set of variables 


requires a knowledge of the second, recourse must bi 
made to some approximate method, which is discusse 
in Section (4). 


(3) SPECIAL CASE WHEN y = | 


This special case which, although it does not appl 
to any real gas, is of interest because it represents 
lower limit to the range of values of y for real gases, is 
an exceptional case to the treatment of the previous 
section. The form (2.1), (2.2) chosen for the com 
patibility conditions becomes singular in this case 
but once a suitable replacement form has been chose 
an analogous treatment may be adopted. This cast 
possesses the further special feature that the char 
acteristic compatibility relations are integrable so that 
exact characteristic parameters are readily availabk 
Unfortunately, because of the nature of the boundary 
conditions, this method does not provide any saving 
of labor in treating shocks of small or medium strength, 
but it may offer a considerable saving if strong shocks 
were to be considered. 

For y = 1 the characteristic form of the equations 
of motion are: on 


dx/dt =u-+c du+t+ (c/p)dp = 0 


on 
dx/dt =u —c du — (c/p)dp = 0 (J.2 
on 


dx/dt =u d&@=0 (3.3 
and the equations of state become 


p = pe*, c=e* (3.4) 
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The shock jump conditions reduce to 
@=0,c =1,u = 6, p = x/(x — 3) 
on a shock moving with a velocity XY given by 
X = sign 6 [(1/2) 6 + V1 + (1/4)6] 


Thus, it follows that an initially homentropic gas with 
= | remains so even after passing through a shock. 


I 
It must be remembered that this has been proved only 


for a perfect gas and that gases with y near unity have 
large molecules and thus show greater deviations from 


perfect gas behavior. 
The characteristic Eqs. (3.1) and (3.2) reduce to 


u + In p = const. = 2aon dx/dt = u +1 (3.5) 


and 
uv — Inp = const. = 28ondx/dt = u — 1 (3.6) 


while the remaining condition (3.3) need be considered 
nolonger. These equations may be treated in a similar 
way to the corresponding ones for steady two-dimen- 
sional supersonic motion,’ and, in this case, the Rie- 
mann function can be determined exactly as a Bessel 
function. Thus, the problem can be reduced to the 
solution of an integral equation, but it is more con- 
venient for present purposes to use an approximate 
method for treating the problem. Thus, as in the pre- 
vious section, a [defined by Eq. (3.5) | and s, which is 
similarly defined to that in the previous case, are used 
as independent variables. Two independent variables 
F and // are introduced by the equations 


dt = Hda+t+ Fds 
dx = H(u — 1)\da + F(u + 1)ds 


where there is now no G appearing since a is a true char- 
acteristic parameter. Moreover, 6 is a function of 
s alone. Thus, a completely analogous, although 
simpler, treatment of the problem can be made yielding 
analogous formulas, and hence there is no need to 


distinguish further the cases y = 1 and y ¥ 1. 


(4) APPROXIMATE SOLUTION FOR THE PRESSURE 
DISTRIBUTION 


To determine the flow field, an iteration process can 
be set up starting from an assumed distribution of 8 
and ® whence the corresponding functions F, G, and 
H can be determined from Eqs. (2.16), (2.17), and 
(2.20). Hence a new distribution of 8 and ® can be 
obtained and the process repeated ad infinitum. If 
the sequence obtained converges then the solution is 
the limit. As a first approximation, the distribution 
postulated in shock-expansion theory is assumed—.e., 
5 = B) and @ = . It can be easily shown that one 
would expect an error in the value of the next approxi- 
mation of the order of &o(%) — x), and this will be 
sufficiently small to neglect in this problem. With 
this postulated distribution, G vanishes and the prob- 
lem reduces essentially to the flow in a homentropic 
simple wave. To reduce the number of sources of 


error, G is eliminated from equations wherever pos- 
sible before the approximate value G = 0 is applied. 
Thus, Eq. (2.27) is to be preferred to Eq. (2.25), es- 
pecially since in the latter the terms, in comparison 
with which G is being neglected, are themselves rather 
small. With this simplification, the differential Eqs. 
(2.16) and (2.17) become 


OF, Oa = Of1,/0s (4.1 


O% fy+1 3 - y¥+1 
(? a+ L *) + = F, = 
Oa \ 2 2 


Ol, (° — ¥ y —1 ) 
a + Dy (4.2 
Os 2 2 


whence, by elimination of //, the first-order linear 
differential equation 
OF, /0a + [n (a — Bo) | Fo = 0 (4.5 


where 


is obtained. This equation has solutions of the form 


Fy = f(s) [(ao — Bo)/(a@ — Bo)]” (4.5 


where the value of f(s) can be made equal to unity by 


suitable choice of definition of s. Thus, s is defined 


by the condition that F = Fy, f = 1 on the shock wave. 
Hence it follows that 
TI, = Fy'(a)s + g(a) (4.6 


where the function g(a) can be determined from the 
boundary conditions on the piston path, which in the 
simplified approximate form become 


Hy = (1/2)T’(a@ + Bo) (4.7 

on 
ds/da = (1/2)T’(a + Bo)/ Fo(a (4.8 
through s = 0, a = a. The piston path can be ob 


tained by direct integration of Eq. (4.8), and g(a@) can 
be deduced from Eqs. (4.6) and (4.7). The approxi 
mate form of Eq. (2.27) which determines the shock 
wave is 
ds di = (j|Fy’(a)] [Fo(a) ]fs + 

) [g(a@)] [Fo(a)}f) [1/w(d)] (4.9 


which can be integrated by the introduction of an inte 


grating factor 


tig n 
(6) = exp | (a — B¥ dé 


As this expression stands, I’ is a function of both 6 and 
do, and this is rather inconvenient. However, in the 
derivation of this equation, 8 has been assumed to have 
the value 6), so that the 6) in the integrand may be re- 
placed by 6 so that T is then independent of 6), which 
simplifies the numerical work considerably. Since 
y ~ y6 for small 6, the above integral is singular at 
6 = 0, but this difficulty is easily removed by writing 


[n/(a — B)] (1/¥) = (y + 1)/(4YY) = 
(2/6) + ¥(6) (4.10) 
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TABLE 2 
y= 1 y = 7/5 y = 5/3 
6 Ir /6° 0/6 I" /6? 0/6 I /6* 0/6 
0) l } l eS ] 3.00 
0.1 1.004 1.03 1.001 3.40 1.001 3.11 
0.2 1.011 $.09 1.002 3.48 1.003 3.21 
0.3 1.025 4.18 1.007 ‘ov 1.006 3.3 
0.4 1.042 1.3 1.011 3.65 1.009 3.47 
0.5 1.067 4.54 1.011 3.72 1.011 3.57 
0.6 1.097 4.79 1.008 3.78 1.012 3.68 
0.7 1.136 5.12 0.999 3.81 1.012 3.79 
0.8 1.182 &.52 0.992 >. 86 1.012 3.88 
0.9 1.235 6.00 0.975 3.85 1.007 3.95 
1.0 1.300 6.60 0.951 3.71 1.001 3.99 


where y, defined by this equation, is integrable near 
6 = 0. Thus, I may be defined by the equation 


(6) = 6° exp | | Yas?" (4.11) 


The equation to the shock then becomes 


l ~ ay (a) 
= 1 Eo" @(8’) do’ (4.12) 
T'(6)J 6 Fo(a) 


where 


0(6) = T'/y (4.13) 


In Table 2 the values of I and © are displayed for 
y = 1, 7/5, and 5/3 for a range of shock strengths. 
Thus, for a given piston path, the shock can be ob- 
tained as a curve in the (a, s)-plane. 

To determine the value 6, of 8 on the piston, use is 
made of the approximate form 

tc) J 574 
By(s) = Bxe(s) 41 + [bd — Px(s)] 

(  4y 


1 &p(s) 
| ays — y(5) ax(s) — [ @(£, s)dé (4.14) 
Y e) a* 


of Eq. (2.28). In the case where entropy variations are 
slow, which numerical work shows holds sufficiently 
well for 6) less than 0.5 approximately, this equation 
can be approximated further to yield 

8,(S) = Bx(s)} 1 + (1 ty) [Bp =a x(s)]} — 

(1/S8y) [ap(s) + ax(s)][®o — Bx(s)] (4.15) 
from which 8, can be calculated immediately. How- 
ever, when this simplification is not available, the dis- 
tribution ®(a, s) has to be determined first. To do 
this, one observes that ® is constant along particle 
paths and hence that 

P(a, s) = P(6) (4.16) 


where 6 is the value of the shock strength at the point 
that lies on the particle path through the point (a, s). 
The differential equation 


ds/da = IIh/ Fo = [Fo’(a)s + g(a)]/Fo(a) 


to the particle path can be integrated to yield 


Ss Sx (6) “ g(é) “o® o£) 
ee ee = a. eee 
Fy(a) Fo [ax (6) | a, [Fy(é) ]? a [F,(§)]? 


The piston path is the particle path with sx = 0, ax = 
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a and hence 


and so the equations to the particle paths become 
s S,(a) + Fo(a)r(6) 1.17 


where 


(6) = } [s — s,(a)]/[Fo(a) ]} $18 


is a known function of 6 for a given motion. Thus. 
Eqs. (4.16), (4.17), and (4.18) enable the distribution 
of entropy to be calculated and the value of 8, to be 


obtained. 


(5) PRESSURE DISTRIBUTION ON A UNIFORMLY 
RETARDED PISTON 


The methods of the previous section have been 
applied to the case of a uniformly retarded piston that 
corresponds to a circular-are airfoil within the accuracy 
of the hypersonic analogy. Graphical methods of 
interpolation and inverse interpolation are sufliciently 
accurate and rapid, and numerical integrations were 
performed using Simpson’s Rule. Solutions were ob- 
tained for each of the three values of y = 1, 7/5, and 
5/3 for the two initial shock strengths 69) = 0.5 and | 
In Table 3, the values of 8, and 8x are tabulated at 
eleven equally spaced times during the piston motion, 
and it is immediately apparent that the variations of 
8, are much smaller than those of 8% (except, of course, 
the case y = 1). In fact, the constancy of 8, is really 
remarkable, and thus, for this problem at least, shock- 
expansion theory is much more accurate than its orig- 
inators could have foreseen. 

In Table 4, a comparison is made between the pres- 
sure distributions on a piston moving with uniform 
retardation as calculated by the three methods; simple 
wave solution, shock-expansion theory, and the method 
of this paper. The value of shock-expansion theory 
over the simple wave solution using the free-stream 
values of 8 and @ is immediately apparent as is the fact 
that the gain in accuracy of the more refined theory 
scarcely justifies the extra calculation in the problem 


considered here. 


(6) ASyMpTOTIC EXPRESSION FOR THE DECAY OF THE 
SHOCK 


In the previous sections, attention has been focused 
on the calculation of the pressure on the piston, and, 
since this is only influenced by the flow nearby (s of 
the order of unity at most), it was not necessary to in- 
vestigate the build-up due to the gradients of 8 and 
®. These terms may become important at large times 
when the shock has become weak. Once the shock has 
become weak, the Friedrich’s Theory will be valid and 
hence will give the correct form for the order of the de- 
cay at infinity, although the coefficients may need 


modification. 
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Calculated Values of B,, 8, 
55 = 0.5 
y= l , = rf 5 = 5/3 
StationT —B, -B, —B, —£, — By —B 
117 Upstream 0 0 2.5 2.5 1.5 1.5 
0 0.0025 0.0025 2.5098 2.5098 1.5118 1.5118 
1 0.0022 0.0022 2 5089 2.5096 1.5108 1.5115 
2 0.0020 0.0020 2.5083 2.5094 1.5100 1.5114 
3 0.0018 0.0018 2.5079 2.5099 1.5094 1.5114 
£18 { 0.0016 0.0016 2.5076 2.5099 1.5090 1.5116 
5 0.0015 0.0015 2.5071 2 5097 1.5086 1.5116 
Thus, 6 0.0014 0.0014 2.5068 2.5097 1.5083 1.5116 
, 7 0.0013 0.0013 2.5065 2.5096 1. 5080 1.5115 
ution 8 0.0012 0.0012 2.5062 2.5095 1.5077 1.5114 
to be 9 0.0011 0.0011 2.5060 2.5094 1.5073 1.5118 
10 0.0011 0.0011 2.5058 2.5093 1.5070 1.5111 
0 = ] 
Y Upstream 0 0) 2.5 2.5 1.5 1.5 
) 0.0188 0.0188 2.568 2 568 1.574 1.574 
] 0.0169 0.0169 2.560 2.570 1.566 1.575 
2 0.0159 0.0159 2.554 2.572 1.559 1.575 
been 3 0.0143 0.0143 2.549 2.572 1.554 1.574 
that | 0.0135 0.0135 2.544 2.571 1.550 1.573 
5 0.0126 0.0126 2.54 2.571 1.547 1.572 
Iracy 6 0.0117 0.0117 2.539 2.571 1.544 1.571 
Ss of 7 0.0112 0.0112 2.537 2.571 1.542 1.570 
8 0.0108 0.0108 2.535 2.570 1.540 1.569 
ntly 9 0.0103 0.0103 2 533 2 570 1.538 1.568 
were 10 0.0100 0.0100 2.5382 2.569 1.537 1. 567 
ob- + The stations of 0 to 10 are at eleven equally spaced times during the uniformly retarded piston motion 
and 
id 1. 
dat TABLE 4 
; Pressure Distribution on Piston 
_— In this table, the pressure on the piston as calculated from the various theories, is tabulated at each of eleven equally spaced instants of 
S ol time during the uniformly retarded motion of the piston 
irse, + =] 
ally : ‘ 
—— ff BN — 
ek. 09 = VO 65 = 1 
rig- Station P(Bao» Po) p(Bo. ® P(Bp, ® pP(Bao, ® p(Bo, ® p(B, ® 
0 1.649 1.641 1.641 2.718 2 618 2.618 
] 1.492 1.485 1.486 2.226 2.144 2.152 
res 2 ] 350) l 343 ] 349 1.822 l 190 ] 768 
3 1.221 1.215 1.217 1.492 1.437 1.450 
rm } 1.105 1.099 1.101 1.221 1.176 1.189 
ple : h C08 0.995 0.997 1.000 0.963 0.975 
6 0.905 0.900 0.902 0.819 0.789 0.800 
10d 7 0.819 0.815 0.817 0.670 0.646 0.655 
TV 8 0.741 0.737 0.739 0.549 0.529 0.537 
- 9 0.670 0.667 0.668 }. 449 0.482 0.440 
im 10 0.607 0.604 0.606 0.368 0.354 0.361 
ict oe 
ry 
0) 1.949 1.9417 1.9417 3.583 >. 463 3.463 
“ml 1 1.714 1.7085 1.7076 2.825 2.746 2.759 
2 1.504 1.49°6 1.4980 2.210 2.160 2.181 
3 1.316 1.3131 1.3134 1.714 1.681 1.702 
4 1.149 1.1468 1.1471 1.316 1.302 1.313 
IE 5 1.000 0.9989 0.9986 1.000 ().997 1.005 
6 0.868 0.8676 0.8674 0.752 0.755 0.761 
7 0.752 0.7515 0.7510 0.558 0.565 0.570 
8 0.648 0.6489 0.6483 0.409 0.418 0.420 
(] 9 0.558 0.5586 0.5579 0.295 0.305 0.307 
4 10 0.478 0.4792 0.4785 0.210 0.219 0.220 
of y = 9/3 
| 
1- 0 2.161 2.1579 2.1579 4.214 $.110 4.110 
qs 1 1.870 1. 8686 1.8670 3.263 3.211 3.220 
2 1.610 1.6112 1.6092 2.488 2.476 2.488 
S 3 1.381 1.3829 1.3834 1.870 1.882 1.882 
s 4 1.178 1.1813 1.1806 1.381 1.408 1.404 
; 5 1.000 1.0040 1.0033 1.000 1.035 1.028 
d 6 0.844 0.8488 0.8480 0.708 0.745 0.738 
4 7 0.708 0.7130 0.7123 0.489 0.525 0.517 
, 8 0.590 0.5954 0.5945 0.328 0.359 0.353 
! 9 0.489 0.4938 0.4929 0.212 0.239 0.233 
0.148 


10 0.402 0.4066 0.4055 0.132 0.153 
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Temporarily, it will be assumed that the piston path 


is known exactly in the (a, s)-plane and that on it // 


Then let 
H=H+h, F=hH+f 


where //) and Fy are the solutions previously obtained. 
Now //, consists of one part that tends to infinity with 
s and another, g(a), which although finite contributes 
importantly to the determination of the shock. Thus, 
although / will be small in the sense that it will have 
an overall scale factor of the order of the maximum en- 


is given. 


(6.1) 


tropy change, it is also likely to possess two terms like 
II, which it will be necessary to split before Eq. (2.27) 
may be integrated. If G is eliminated between Eqs. 
(2.16), (2.17), and (2.20), the resulting equations, in 


terms of f and / become 


Of Oa = Oh Os (6.2) 
and 
Oo|.3- 7 

Fy (3 7 Bo) bh 
Oa 2 
(TI, + h) (Of Os) + | 

i + Le? a + 
{ (u +c) [2y(7 — l)c — (08 Oa) | 2 


3 - : 3) 
8) - 2]? ¥ ! Ho(8 — &) +h X 


2 
7-4 : a. a 
( 2 at 9 a) (6.:3) 


These equations have to be solved subject to the 


boundary conditions. 
h=Qons = Sp(a@) 
and f = 0 on the shock wave, which is given by 


ds/dé = [UIy + h)/ Fy] }1/[W(6)]} (6.4) 


(with s = 0, 6 = 69). Since the isentropic orders of 
magnitude can be expected to apply for large s and, 
therefore 6 = O(s~'’), Eq. (6.4) implies that # can be 
O(s) at most. Thus the asymptotic form of (6.5 Eq.) 
becomes 


(Of/Oa) (a — B) = Bo) + O(®mars % 


(OFy/Oa) (8 — 
from which it can easily be deduced that 

h = [(B — Bo)/(a@ — B)] (OF)/O0a)s + R(a, s) (6.5) 
where k is a bounded function of order ®). Thus, the 
Eq. (6.4) to the shock can be written in the form 


ds 

dé 

g(a) + Fy’(a)} 1 + [(8 — B)/(a — B)|t s + Rk(a, s) 
Fo(a) (6) 


whence it follows that 


(ds/d6) + [(y + 1)/(4 Vy) ]s = 
[g(a) + Ra, s)]/[Fo(a)y(6) | 


(6.6) 


and hence that the exact equation to the shock is given 


by 
io my gtk 
T'(5)J om Fy 


(6.7) 


0(6) dé 
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It is to be observed that this solution is consistep 
with the postulated Friedrichs’ form assumed in « 
riving it. If the appropriate expansions are now mac 
s = [K(6o)/6?] [1 + 0(67)] [1 + O(@ 


where 
K i) 
To convert this asymptotic expansion to one involy 
ing the physically measurable variable ¢ rather thay 


s, Eq. (2.14) is 


integrated along the shock and hence one obtains the 


{= | Fila [ + W(6 AL 
J0 ao 


Thus, one can obtain the asymptotic formula 


the mathematically significant variable 


equation 


t = Fo(a:)s{1 — Ak ’’s 1 () 
where 
= ae [/. 1 a3 
A = V\ 1 a oa l } || / — i ct r 7 T | 


which implies that A is known exactly, even allowing for 


the effects of entropy gradients. Thus, the decay of 
the shock at large times is given by 
§= ft | Fo( an) K (69) fi — (1, 2)A;¢t- 

[ Fy(a1)K (6) ]} + O(¢-')] [1 + O(%))] (6.8 


which is to be compared with the Friedrichs’ solution 
6= [4 (y+ 1] [VA et —-— (4 t) + 0 )] 

F,K (69) 1s 

evaluated under the assumption that 6, is small, then 


when the shock is everywhere weak. If 
Eq. (6.8) reduces to the Friedrichs’ expression. 


(7) EXTENSION OF RESULTS TO PLANE FLOW 


The usefulness of shock-expansion theory in dealing 
with comparatively strong shocks in one-dimensional 
unsteady flow has been amply illustrated by the results 
of Section (5). The important question remaining to 
be answered is whether a similar result holds for two- 
dimensional supersonic flow, and, if so, for what range 
of shock strengths and Mach Numbers. Certainly 
from Goldsworthy’s work on the “hypersonic analogy” 
the result must be true for sufficiently large .1/ and cor- 
respondingly small leading-edge deflection. However, 
some doubt must remain as to whether this will cover 
completely the range of practical importance. 

The method developed in this paper may be extended 
to the two-dimensional case, but the equations are not 
so tractable and there is a further parameter, the 
Mach Number, involved. However before such an 
analysis is undertaken it is pertinent to inquire under 
what conditions the tendency of errors to annul occurs. 
If the appropriate form of the characteristic equations 
is examined, then it can be shown that for sufficiently 
small stream deflections the errors tend to cancel if the 


(Continued on page 720) 
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On Certain Second-Order Effects in the 


Limit Design of Frames 


E. T 


Brown 


ABSTRACT 


The second-order effects due to the changes in geometry are 
nvestigated for the frames composed of plastic-rigid material 
In the case of ideally-plastic material, a criterion is developed to 
determine whether the quasi-static distortion of the frame pro- 
eeds under increasing or decreasing loads. The analysis also 
shows that when the compressive axial forces are negligible even 
. slight amount of strain-hardening considerably increases the 


iccessary for quasi-static deformations 


load 1 


(1) INTRODUCTION 


HE PLASTIC METHODS OF ANALYSIS and design for 
| po for beams and frames" * are mainly concerned 
with the load at which the structure could become a 
mechanism following the formation of plastic hinges 
at a sufficient number of sections. 

If the changes in geometry caused by the subsequent 
deformation is neglected, then it can be shown that in 
the case of frames of ductile material the deformation 
How- 


ever if this change in geometry is taken into account, 


of the structure proceeds under constant loads. 


then continuing deformation under constant external 
loads is possible only in exceptional cases.* 4 Ac- 
cordingly as a rule, quasi-static flow then requires 
If the 


applied loads are increased monotonically, as is cus- 


either increasing or decreasing external loads. 


tomary, the first case corresponds to a stable process, 
whereas the second case corresponds to sudden col- 
lapse. Since the present first-order theory of limit 
analysis does not furnish a criterion that would allow 
one to distinguish between these two cases, it seems 
worthwhile to develop a second-order theory that 
should furnish a clear criterion for collapse. 

If the frame is composed of a rigid strain-hardening 
material in which Young’s modulus is infinitely large, 
theorems of limit analysis define a load at which de- 
formation first becomes possible. In_ this the 
load increase to enforce further deformation is governed 


case 


by the rate of strain-hardening and also by the changes 


of shape. Because of strain-hardening, the plastic 
deformation in the frame is no longer concentrated at 
isolated sections, and, instead, this deformation is 
spread over certain finite segments of the frame. 

The purpose of the present paper is to determine the 


rate at which the loads need to be increased to enforce 
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a quasi-static motion of framed structures made of 
either ideally-plastic rigid or rigid-plastic strain 
hardening material. 

It will be shown that in the case of ideally-plastic 
frames an energy criterion can be established to dis 
tinguish between the two cases mentioned above. In 
the case of strain-hardening materials, it is found that, 
if the influence of the axial forces is neglected, then the 
subsequent deformation can proceed only with increas- 
ing load when the rate of strain-hardening is not very 
small. 

Asa note of caution, it must be added that the second 
order analysis developed herewith must be supple 
mented by a study of elasto-plastic stability whenever 
compressive axial forces of significant intensity exist in 


the structure. 


(2) FRAMES OF IDEALLY PLASTIC KIGID MATERIAL 


Let us consider a frame composed of ideally-plastic 
rigid material, and suppose that the loads on this struc 
ture are steadily increased in fixed proportion to each 
other. A critical load is finally reached at which the 
fully plastic moment .1/p is developed at a sufficient 
number of sections, and the structure becomes a mecha 
nism. The calculation of this critical load and the corre 
sponding mechanism is fairly simple in comparison 
with elastic analyses for highly redundant structures. 
The subsequent quasi-static motion of this mechanism 
and the rate of change of the loads to achieve this mo- 
tion constitute the scope of the present investigation. 

In order to fix ideas we shall first consider a very 
simple model. 

The portal frame shown in Fig. 1 is composed of mem 
bers having the fully plastic moment 1. It 1s easily 
shown that when P reaches the critical value 4.1/) 
plastic hinges develop at the corners A, B, C, D, and 
the frame becomes a mechanism with a single degree 
of freedom. The subsequent motion of this mechanism 
is governed by the basic assumption which states that, 
if at any one cross section the moment retains the 
value .\/,, then this cross section acts like a hinge (see 
Fig. 2 

This suggests that the quasi-static deformation of the 
considered mechanism is only possible if the bending 
moment remains constant, and equal to W/o, at A, B, C, 
and D. 

The resulting four conditions, together with the three 
equilibrium equations, yield the following equation for 
the instantaneous value of P as a function of @ (see 


Fig. 1 
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P = (4M,/l) (1 + 4) = 
(4My/1) (V 2/2) cosec [6 + (7/4)] (1) 


and therefore for small angles of rotation occurring im- 
mediately after the yield point has been reached, 


It follows that quasi-static motion can only be main- 
tained with decreasing loads, and so P = 4./)// is the 
load at which real collapse occurs. 

However, if the loads acting had the reverse direc- 
tions as depicted in Fig. 3, then the same analysis 
shows that the quasi-static flow would require an in- 
crease in the magnitude of loading, the instantaneous 
load now being given by 


P = (4M)/1) (V 2/2) sec [6 + (2/3)] (2) 


Load deformation curves for both cases are shown in 
Fig. 4. The influence of the axial forces limits the 
validity of Eqs. (1) and (2) to the range where @ is small. 
The above analysis can easily be extended to more 
complicated cases. However, in general this will in- 
volve the rather lengthy elimination of internal forces 
from the equations. In the next section it will be 
shown that an alternative and simpler procedure is 


available. 


(3) ENERGY CONSIDERATIONS 


We shall now establish an energy criterion for dis- 
tinguishing between our two collapse cases. In the 
example of the previous section let us compare the 
energy dissipation I]! with the work V done by the 
original loads after a deformation corresponding to a 


rotation @. If for small rotations 


W(@) > V(@) 
we may conclude that quasi-static deformation re- 
quires additional work and hence increasing loads. On 
the other hand, if 

W(é) < V(@) 
then the loads must decrease. Since 


W(6) ~ 4110 / 


we 


(0) ~ 4, [0 + (@2/2)]\ 


the positive sign being taken in the first case (Fig. 1) 


and the negative sign in the second case (Fig. 3) 

Vo) — W(@) ~ + 2M (4) 
which is in agreement with the conclusions based upon 
the statical considerations mentioned earlier. This 
procedure can be generalized without difficulty to more 


complicated cases. 


(4) FRAMES OF RIGID STRAIN-HARDENING MATERIALS 


For the purpose of the present investigation, it is 
sufficient to consider a frame composed of material 
having the idealized uni-axial stress-strain diagram 


1955 


OCTOBER, 


shown in Fig. 5. For continuous loading in tensjoy 
or compression, we have then the following stress-strain 


relation; 


oO = | € + do 5 


where all symbols are defined in Fig. 5. 

If the cross section of the beam is doubly symmetri: 
and if the bending moment lies along an axis of sym 
metry, then the classical assumption that the cross 
sections remain plane during bending yields the fo] 
lowing relationship between the bending moment ]/ 
and the radius of curvature p of the neutral axis 


(M — M)/(E,1) = 1/p when | M) > | MM) ) 


1/p = Owhen | M| < |M \ 


where F, is the plastic tangent modulus of the ma 
terial, F is the cross-sectional moment of inertia, and 
My is the initial yield moment. 

Now, let us investigate the state of affairs in the 
neighborhood of a yield hinge H following the in- 
cipient motion (see Fig. 6). Suppose that during the 
subsequent motion the bending moment at the cross 
sections in the neighborhood of H increases. 

According to Eq. (6) the neutral axis of the beam will 
be curved in the region AB, and this curved part of the 
beam will generally change its size and shape during the 
deformation. 

Since at the end points of the curve AB 


l/p = 0 


the bending moment computed on the straight line 
extension AH’ and H’B of the members | and 2 can 
reasonably be expected to be sufficiently close to the 
actual bending moment. This suggests that the initial 
deformation of the frame is closely approximated 
by the deformation of a frame having ideal yield hinges 
However, the simple yield condition for ideally-plastic 
frames should be replaced by an equation expressing 
the relationship between the angle of relative rotation 
6,4, Of two neighboring rigid portions and the bending 
moment distribution. 


From Eq. (6) 


aaa f- 7 [ M—- M, F 
7 ae p Ja E,1 - 


and therefore since the bending moments are computed 
along straight lines AH’ and H’B, it follows that 


M mas 


| 
Elan = Y g (Mo — My)d MV + 
le 


l *Mna 
= | (Ml — Md 


VW 


or 


E,l0an = [(1/Vi1) + (1/V2)] (Minar — Mo)? (i 


where V; and Il’, are the absolute values of the shear 
forces (Fig. 6). 

Equation (7) is the condition for strain hardening 
which replaces Ming: = Wy for perfectly plastic frames. 





SECOND-ORDER EFFECTS IN THE LIMIT DESIGN OF FRA MES 683 





-nsion 
Strain 


p C M 4M, /2! 


netri P BY — 4 
sym | | 
/ 
: 
y 

















CTOSS 
a fol- 48 
nt J/ | 























* Ma- 


, and Fig. i Fig. 2 Fig. 3 


n the 2 M 4 (2) 
- [F) 


Cross 


n will 1.0 
f the 
ig the 


— 
— 


Cia -steas iny. eacias ee 


- line 


2 can 


o the ] l ! | | —_—_»> be 
nitial 0 10° 20° 30° 40° 50° 60° 
nated 


nges. Fig. 4 

lastic oC 

Ssing ‘ | P 
H 


ation 
iding ear 
| 
Do | ra I 1 I TT 
wont TULL |M: FLEE 
aii M 
a ia Ss H’ . 


eee ee ee TT 






























































shear Fig. 5 Fig. 6 


ning 
mes. 



































684 JOURNAL OF THE AERONAUTICAL SCIENCES OCTOBER, 1955 
a and H becomes a yield hinge. If the angle of rotatioy 
Load = P/P, ; : aera i e 
ars s y 5 oO 
| of one of the bars is denoted by 6, the condition (7 
H8 applied at H gives 
2F,10 = [Pl cos (a — 0) — My]?/[P cos (a — 6)] 
1.2 : 
or 
oe 
PP, = [cos a/cos (a — 6)] X 
1.0 1 + 6+ [O(2 + 6)]*} (9 
deally p 
| j €70 
OF °Stic where 
(E,1)/(Mil) = « 
The corresponding load deformation curves are shown ™ 
in Fig. 7 for a typical shape of frame and for various Num 
degrees of strain-hardening. It is clear that even a zles | 
slight amount of strain-hardening will considerably ant 
° ° . ° e 7 tur 
SL increase the load necessary for quasi-static deforma- 3 i 
tions of the structure. ss 
For example, if the frame cross section is rectangular for « 
tinuo 
2 e = (1/3) (E/o») (h/l) the ke 
: . ' for tl 
where / is the depth of the bar. For the relatively ont 
small degree of strain-hardening E,/o) = 1.5, and neces 
fori /k = 10, ¢ = 0:05. istic 
Whereas in the case of ideally-plastic material P = lengt! 
fe) 4 4 4 4 > > : + - ‘ flectic 
.05 A AS 2 Py is a real collapse load, a slight amount of strain rt 
. - . . . . 0; 1€ 
@(deformation) hardening of the material necessitates increasing loads feels 
with increasing deformations (Fig. 7). 
Fic. 7 
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Design of Two-Dimensional 
Continuous-Curvature Supersonic Nozzles 


JAMES C. SIVELLSt 
Sverdrup & Parcel, Inc. 


ABSTRACT 


In order to cover continuously a rather wide range of Mach 
Numbers, several supersonic wind tunnels have flexible-plate noz- 
zles incorporated in their design. The elastic curve of a flexible 
plate supported at discrete points must have continuous curva 
It is, therefore, desirable that the curvature of the aerody- 


ture 
the flexible plate 


namic contour also be continuous to enable 
to conform to the desired contour. 
for designing two-dimensional supersonic nozzles having con- 


are made for determining 


A method is presented herein 


tinuous curvature. Approximations 
the length of the nozzle and for correcting the perfect-fluid profile 
for the effect of the boundary-layer growth. Computing forms 
are included for making the successive approximations that are 
necessary to establish the angles of the inflection and character- 
istic points that will satisfy the desired requirements of nozzle 
length, throat radius, and continuous third derivative at the in 
flection point. After these angles are established, the coordinates 
of the perfect-fluid profile, downstream of the inflection point, can 


ye obtained by the characteristics method. 


SYMBOLS 


1, B <&,. % 
O, P,Q, R, 
S, S’, T, U = points on the characteristic network of the 
nozzle as defined in Fig. 2 
= flux fraction [see Eq. (19) 
= area 

= = throat area 

b = flux fraction [see Eq. (27) 

G = curvature 

D = distance from throat to the end of the flexible 
plate 

H = physical half-height of test section 

H = ratio of boundary-layer displacement thick- 
ness to momentum thickness [see Eq. (41a) 

h = one-half test section height before boundary- 
layer correction 

R = an integer 

K = anangle defined by Eq. (7) or (8) 

length of nozzle measured from the throat to 
the first point of parallel supersonic flow 

L = overall length of the nozzle measured from 
the throat to the point of cancellation of the 
last expansion wave 

Mf = Mach Number 

m = an integer 

\ = boundary-layer velocity— profile parameter 

n = a number, defined by Eq. (34) 

R = Reynolds Number [see Eq. (4la 
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R* = radius of curvature at the throat 

R’* = radius of curvature at the throat corrected for the 
effects of the throat boundary layer 

? = an integer 

= coordinate along contour APQ, positive to right 
(see Fig 2) 
coordinate along AO, positive downward 

v = coordinate along centerline of nozzle 

} = coordinate of the initial region corrected for bound- 
ary layer 

y = coordinate normal to the centerline of the nozzle 

y* = coordinate of throat normal to the centerline of the 
nozzle 

6 = boundary-layer displacement thickness 

7] = flow inclination angle 

Oy, = inclination of flow along contour 

A; = inclination of flow crossing AO 

m = Mach angle 

v = Prandtl-Meyer angle 

¢ = slope angle of characteristic lines defined by Eqs. (23 


and (29) 


Subscripts denote point on characteristics network 2 


See Fig 


INTRODUCTION 


— DESIGN OF TWO-DIMENSIONAL supersonic nozzles 
by the method of characteristics has been treated 
in numerous papers and books. The general supersonic 
nozzle may be divided into two regions: the initial 
region (from the throat to the inflection point) where the 
nozzle surface is gradually curved away from the 
centerline, and the terminal region (from the inflection 
point to the test section) where the nozzle surface is 
curved back until it is parallel to the centerline. In 
order that the test shall have shock-free, 
parallel flow, the expansion waves that are created and 
reflected in the initial region must be canceled exactly 
by the compression waves that would be created in the 
terminal region if they did not cancel the expansion 


section 


waves. 

In the general design of a nozzle by the characteristics 
method, the contour of the initial section has often been 
arbitrarily chosen. The contour then deter 
mines the shape of the terminal curve. 
method is to choose a transverse Mach 
stream angle distribution at the inflection point and 


initial 
An alternate 
Number and 


then use the characteristics method to determine both 
the upstream and downstream contours. In either 
method, the determination of the contours is often done 
graphically, although the accuracy can be improved 
by using analytic geometry methods to compute the 
locations of the the characteristics. 


The accuracy of the final curve increases as the strength 


intersections of 
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Fic. 1. Second derivative of a Foelsch nozzle 


of the characteristics used decreases, which, in turn, in- 
creases the number of intersections to be computed. 

By assuming that the flow at the inflection point is 
radial—i.e., that it originated at a single source 
Foelsch! showed that the terminal curve can be com- 
puted analytically. The initial curve must still be de- 
termined graphically or by analytic geometry, although 
in some cases a simple analytic curve (a cubic in ref- 
erence 2) has been used to approximate the initial 
curve. The use of such an analytic curve obviously 
does not ensure radial flow at the inflection point. 

In the design of the majority of supersonic nozzles in 
use at present, the characteristics throughout the 
length of the nozzle are all of equal strength. In the 
initial region, the creation of the expansion waves de- 
mands that the curvature (of the upper contour) have 
a finite positive value. Similarly, in the terminal region, 
the cancellation of the expansion waves demands that 
the curvature have a finite negative value. If a pure 
reflection occurs at any point, the curvature is zero. 
This means that the curvature is discontinuous, par- 
ticularly at the inflection point and at the test section, 
These discontinuities can be readily calculated for 
Foelsch’s analytical curve (Fig. 1), but they are also 
present in any design wherein all the characteristics 
have equal strength. Such discontinuities can be 
machined into solid block nozzles but cannot be matched 
by the elastic curve of a flexible-plate type of nozzle. 
The inability of the elastic curve to match the theoreti- 
cal aerodynamic curve produces undesirable variations 
in Mach Number and stream angle in the test section. 

Since the curvature of the contour must pass through 
zero at the inflection point (by definition), an aerody- 
namic curve must be found which has continuous curva- 
ture through this point if it is to be matched by an elastic 
curve with any desired accuracy. To improve the 
accuracy of the matching, it is desirable that even the 
third derivative be continuous through the inflection 
point since it is impractical to maintain the inflection 
point at a given support point for any reasonable range 
of test section Mach Numbers for which the flexible plate 
nozzle may be designed. Norman Riise of the Jet 
Propulsion Laboratory, California Institute of Tech- 
nology, Pasadena, Calif., working with Dr. Allen E. 
Puckett, developed a method of designing a continuous 
curvature supersonic nozzle.* This method utilizes a 
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transition region following the inflection point, in which 
region the characteristics are partially canceled and 
partially reflected in a manner to make the third deriyg 
tive of the contour continuous in this region. Following 
this region is the terminal region in which the charac 
teristics are exactly canceled including those which wer 
partially reflected in the transition region. The method 
is still that of characteristics but is one that is readih 
adaptable to automatic computing equipment. 

In any method of obtaining contours for a supersoni 
nozzle, it is desirable to be able to determine the overal] 
length of the nozzle before obtaining the ordinates 
Also, the right combination of variables is required to 
insure that the third derivative of the transition curve 
is equal to the third derivative of the initial curve at the 
inflection point. This report presents a method of ob- 
taining this right combination of variables together with 
an approximate value for the length of the nozzle. Inas- 
much as the ordinates of a flexible-plate nozzle can be 
corrected to a certain amount during calibration tests, 
the boundary layer growth does not need to be ac- 
curately calculated. Some simplifying assumptions are 
given herein for taking into account the boundary-layer 
growth in obtaining the actual ordinates from the free- 
stream ordinates resulting from the characteristics 


solution. 


METHOD OF OBTAINING CONTINUOUS CURVATURE 


The method of calculation derived by Riise’ is 
briefly outlined below. Referring to Fig. 2, the point A 
is the inflection point. The assumption is made that the 
flow is radial, so that the are AO, of a circle around the 
apparent source, is an equipotential line along which the 
Mach Number is constant. 

Any point along the contour APQ, s,,, is defined by a 
corresponding point /; lying on the same characteristic 
If all characteristics from point A to and including the 
point s, were exactly canceled at the wall, then 


0. = (0, + 6,)/2 


u 


and 
d6,, = (1/2)d@, (Z 
These relations are true even though the flow is not 
assumed to be radial, so long as the Mach Number at 
t, is the same as that at A, since 9, — 0, = 0, — &. 
The curvature of the contour just to the right of point 
A is 
Cs = (6B;/dse) a | 
(1 /2)(d0; /dt,)4(dt,/ds,)4 (3 


Il 


(1/2)(d0,/dt;)4(1/V M4? — 1) 


The value d6, dt, has the magnitude of the curvature of 
the arc AO but is negative. Even though the flow were 
not assumed to be radial, d6,/dt,; would be finite and 
negative and therefore the curvature C, would be finite 


and negative. 
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If any of the characteristics between points A and s, 
were not exactly canceled but were partially reflected, 
§, — 0. < O, — 6, or 


6, yo (O64 + 6;) 2? 


The rate at which 6, and therefore the curvature, varies 
can be specified. Riise and Puckett chose the relation 


d6,, = |(0, — 0)/K]d@, (4) 


where A is a constant, but many other relations could 
have been chosen that would accomplish the same result 
of making (d@,.)4 and therefore C4 equal to zero when 
#, = %. Integrating the Riise and Pucket relation 
gives 


0. = O04 — [(01 — 0,)?/2K] (5) 


The value of A can be determined by selecting a point 
P beyond which it is desired to completely cancel the 
incoming characteristics. The characteristic to be can- 
celed at P intersects the are AO where the flow inclina- 
tion 1S 9; p. 
Therefore 


(d0,.)p = [(04 — O@p)/K ](d0,)p = (1 °2)(d0;)p (6) 

or 
6, — Op = K/2 (7) 

also 
6, —Op = K'S (8) 


Thus, A, 6;,, and 6p are interrelated, and any one of the 
three can be considered to be the parameter defining 
the rate of change in the curvature. The point P has 
been defined as the characteristic point’ since the 
characteristic from P to T (Fig. 2) defines the beginning 
of the region of uniform flow at the point T. This 
characteristic PT and also TU have zero strength and 
therefore d.\J/dx = 0 along characteristic TU, and the 
curvature of the contour at point U is likewise equal to 
zero. The characteristic point must always lie down- 
stream of the inflection point. In order for the charac- 
teristic AR to have zero strength, the contour must have 
additional inflection points upstream of A (see reference 
dD). 

The value of the Prandtl-Meyer angle at the test sec- 
tion, designated vr, is described by the test section Mach 
Number. The value at the characteristic point is 


vp = vr — Op / 
(9) 
= yp — 0, + (K/8) \ 
The value at the inflection point is 
vis = v~p — (Op — Hp) | 
» (10) 


= yr — 0, — (K/4) \ 
he design of a supersonic nozzle by this method re- 
quires specification of both parameters 6, and A (or 
%,,0r 6p). Inasmuch as 


V4 > 64 (11) 








bo 2 4 
}. L | 
Or Ox Va + Vr-Ga- tan Pas = a’ tan din + (>a ) ran (pt Ge) 
Gr: Qu-"Ye Vr ° Vy-% van Ps1 = b TAR (ps-6q +(I- b’) ram Ad ’ 


Os. }s «)1-% 


Fic. 2. Characteristic diagram for a continuous curvature 
nozzle 


the limiting value of @4 is given by the relation 
0, < (vr/2) — (K/8) (12 


The parameter A can have any value from zero (Foelsch 
contour) to the value of 40, (when @;, = —6,). From 
a practical standpoint, it is believed desirable to keep 


the value of K above 26,. 
APPROXIMATE LENGTH OF NOZZLE 


Terminal Length 

The following method of estimating the terminal 
length of the nozzle gives much more accurate results 
than that given in reference 3. Referring to Fig. 2, the 
flow has been assumed to be radial, normal to the arc 
AO. The characteristic AR defines the downstream 
limit of this radial flow. 

The following relations hold: 


arc AO/y* = A4/A* = (2/180) (6474/y*) 
(6, in degrees) (13) 
h/y* = Ar/A* (14) 
(7/180)(@4r7r/y*) = An A* (seereference 1) (15) 


The Mach Numbers at A and R are obtained from the 


values of v4 and vp where 
va = vr — 6, — (K/4) (16) 
VR = VT (K 4) (17 
Therefore 


length OR = rr — Pr 
(180/7)(h 64) } 


II 


(Ap/A*) — (A4/A*)]4 
(Ar A*)! (18) 





This relation is exact and does not involve any ap- 
proximations. 

The determination of the length RT involves some 
approximations. The Mach Number at S can be ac 
curately determined from the relation vs = vr — (A/S) 
and the stream angle at S is equal to A/S. The point 
S is the characteristic point on the streamline for which 
point C is the inflection point. The characteristics 
crossing CS have the same strength as those hitting and 
reflecting from the contour AP. Therefore 
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0c — Os = 04 — Op = K/8 


and 06¢ = K/4. 
16,;, the points C and A merge, the points B and R 
merge, and the points, P, Q, and S merge.) The flux 
crossing the characteristic RS at Mach Numbers varying 
from J/7p to Ms is that which passes through the angle 
K 4 at M4. It is assumed that this flux can be divided 
into two parts: one part cross the the line RG, drawn 
tangent to RS at R, at a Mach Number of JA/p and the 
other part crossing the line GS, drawn tangent to RS at 
S, ata Mach Number of l/s. The point G is shown in 
Fig. 2 as the intersection of these tangents. Also shown 
is the point S’ which is the projection of the point S on 
the axis. Let a be the fraction of the flux at the Mach 


Number of |/p. Then 


RG sin MR (Ap A) zs lamr4(K 1) ] (1SOA 4 ma) (19) 


GS sin us/(As/A*)=(1 — a) [mr4(K/4)]/(180A 4 /A*) 


(20) 
Furthermore 
S/S = ys = RG sin wr + GS sin [us + (K/8)] 
wra(K }) | Ar 
= a 1 
1S0A 4 /A* a” (21 
: As(_K org =) 
( a) A* cos 8 + vV Ms? — sin Q 
RS’ = RG cos pwr + GS cos [us + (AK/S)] 
mrs(K 4) | Ap “fn ee ” 
= a y . et —<¢Z) 
1804, A*|° A* = (22 


As K TY I 
V Ms? — 1 cos — sin 
An S 8/ | 


This treatment is analogous to defining the slope of the 
straight line connecting R with S by the relation 


tan drs = (S’/S/RS’) = a’ tan wr + (1 — a’) X 


tan [us + (A'8)] (23) 
where a’ is related to a by 
a (1-—a@) 
a 7 a 
As/A*[V Ms? — 1 cos (AK 8) — sin (K/S)] (iia 


Ar/A*V Mr? — 1 


If the characteristic RS were a cubic, a’ would have a 
value of |; since the curvature of the characteristic is 
zero at the point S$. A comparison of several nozzle 
contours obtained by the method of characteristics in- 
dicates that a very good approximation for the lengths 


S’S and RS’ is obtained when 


(Ap/A*)V Mp? — 1 


(As/A*)[V ‘Ms? — 1 cos (K/8) — sin (K 8) ] 
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(When A has its maximum value of 
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1 — a’ 
: = 7 <6 
a 
/ , ; ; 
(As/A*[V Ms? — 1 cos (K/8) — sin (K/S8)] 
26) 
(Ar/A*)V M,? — 1 


An extension of the above treatment to the charac. 
teristic ST leads to the relations 
rrs(K 4) 


S'S = x 
180A , /A* 


| As ( K ) : =) Ar] 
6 — |cos — VM.2 — 1 sin + (1—5 
A* 8 8 A* 


S/T = mra(k }) , 4s " - K R 
~ 180A, /A*  A* Ms’ — 1 cos Qt sin) 


A* 
and 
tan gsr ~ b’ tan [us — (K/8)] + (1 — 5’) tan yy 
(29 
where 
1 — 0’ (1 — bd) x 
yo b 
(Az A*)\* Mr? — |] 
(30 
(As 'A*)[V Ms? — 1 cos (K/8) + sin (K/8)] 
and 
p= 
a(Ar/A*) + (1 — a)(Asg/A*) [cos (K/8) + 
V Ms? — 1 sin (K/8)] — (Ar/A* 
(As A*) [cos (K/8) — V Ms? — 1 sin (K/8)} —- 


inasmuch as both Eqs. (21) and (27) for S’S must result 
in the same value. 
Combining the above relations, the length Ax is ob- 
tained by the relation 
A ISO jAr A, 
os = ii ; — , COS | ten 
70;(Ar A*) (A* Pag 
(A 4) Ar 
a 
180 | A* 


Pd 
= 


As ( K ; K 
a ~ | V VWs? — 1 cos sin ) 
A* : 
A s A 
b 1+ V Ms? — 1 cos + sin 
ij i! 9 
(l —O 1* V AM; = ac 2 
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DESIGN OF 


x = Avr+ iV M7? — 1 


where the approximation involved is in the value of a 
which determines the length RT. This length RT be- 
comes an appreciable part of the terminal length when 
K @,islarge. The range in values of a has been found to 
be from ! to */s, depending upon Mach Number, 


K @,, and 64. Increasing the value of A 64, or decreas- 


ing the value of 6, for a given test section Mach Num- 


ber, increases the overall length of the nozzle. 


INITIAL LENGTH 


Inasmuch as the use of a flexible-plate type of super- 
sonic nozzle allows some correction of the ordinates to be 
made during the calibration period, the design of the 
upstream portion can be simplified by using an analytic 
curve from the throat to the inflection point. The use 
of such a curve has been found to give nearly radial flow 
at the inflection point and circumvents the need for de- 
termining accurately the boundary layer in this region 
and the shape of the sonic line at the throat. If the de- 
sign must encompass a wide range of Mach Numbers, 
the length required for the higher Mach Numbers 
may be too long for the lower Mach Numbers even 
though the parameter A is allowed to reach its maxi- 
mum value of 40,. Under these circumstances, either 
the throat position must be moved downstream for the 
lower Mach Numbers or the initial region must be 
lengthened. An analytic curve that allows a long initial 
region 1S 


y = y* + Ax? + Bx" 34 


The conditions that this curve must satisfy are: 
|) The ordinate must be equal to y* when x = 0. 


2) The slope must be zero when x = 0. 
3) The ordinate at the inflection point must be equal 
tor, sin 6,4. 
1) The slope at the inflection point must be equal to 


tan 64. 

5) The curvature at the inflection point must be 
equal to zero. 
When these conditions are satisfied the expression for 


the initial curve becomes 


(mn — 1) tan 6, : l tan 6, 


2in — 2) Ny n(n — 2) x4. 


or 


(n — 1)x4 tan 6, (* 2 x, tan @, [x \" 
agi 7 
. 2(n — 2) x4 nin — 2) \xy 


36) 


Where x4, the length of the initial curve projection on the 
axis of symmetry, is given by the relation 


2n h}(A4/A*) [180 sin 64 /(764)] — 1! 


V4 = 
(A7/A*) tan 0,4 


(n + 1) 
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Nonintegral values of m may be used with these equa 
tions. It can be seen that when 7 = 3 the equations re 
duce to the cubic used by Crown.” In order to keep the 
overall length of the nozzle to a minimum, consistent 
with other requirements, it is suggested that the cubic 
be used when the value of A is less than its maximum 
value. 

One consideration that must be taken in the struc 
tural design of a flexible-plate nozzle is to maintain the 
radius of curvature at the throat above a minimum 
value. The radius of curvature obtained from the gen 
eral analytic curve is 


R* = [(n — 2)x4]/[(n — 1) tan 64] 38 


In general, the throat radius of curvature will be a 
minimum for the maximum Mach Number for which 
the nozzle is designed. The minimum value acceptable 
structurally is usually the determining factor in the 
selection of the parameters 6, and A, on which the over 


all length of the nozzle depends. 


CORRECTION TO CONTOUR DUE TO BOUNDARY LAYER 


An examination of the growth of the boundary layer, 
calculated for a number of supersonic nozzles, indicated 
that the growth is approximately linear from the inflec 
tion point to the test section. Inasmuch as an analytic 
curve was arbitrarily chosen for the initial region, an 
analytic curve describing the growth of the boundary 
layer in this region can be chosen just as arbitrarily, if 
such a curve satisfies the conditions: (1) the displace 
ment thickness at the inflection point has the right 
value and (2) the rate of change of displacement thick 
The addi 


) 


tional conditions are desirable mathematically: (3 


ness is continuous at the inflection point. 


the curvature of the displacement thickness curve at 
the inflection point is zero and (4) the rate of change of 
displacement thickness at the throat is zero, This 
fourth condition is necessary to avoid moving the posi 
tion of the geometric throat. The usual assumption is 
made that the boundary layer has zero thickness at the 
throat. Actually, the boundary layer may not be zero 
at the throat, and the aerodynamic throat may be 
downstream of the geometric throat. The displacement 
thickness at the inflection point, will, however, be ap 
proximately right. 

The above conditions lead to the following analytic 


relations for the displacement thickness. For 0 <x <x, 


6*/64* = (3/2)(x/x4)? — (1/2) (x/x4 39 
For vA <x < D 
6*/ip* = (3x — x4)/(83D — xa 10) 
where 
6,*/bp* = 2x4/(383D — x4) (41) 
It should be noted that the distance D is used 
instead of L because the end of the flexible plate is 


where the physical half-height is known. The distance 
L may be smaller or larger than D, depending upon the 
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TABLE | 
Terminal Length of Nozzle 
Mr = 3.0, vr = 49.7573, Ar/A* = 4.23457, /Mr? -— 1 = 
2.82843243, (Ar/A*)/ Mr? — 1 = 11.9771951 

[a] & 7» 7057 
[el 7 : 
| 3) K/4 = (1) xf2)/4 7+ 7057 
[@! w= v-O-@ _ Boe 3459 
| (5) M, FROM TABLES @ @) 2. 302518 
[| As/a* FROM TABLES @ ©) 2.198143 
[@| COS &= cos WO 0» 990969862 
[ ~e| Ye = Vr - | 42.0516 
[@|_ Mag FROM TABLES @ 8) | 2.628476 
LG] Vk -1 = AQ | 2e4308200 
(42)! Ag /ax FROM TABLES @ | 2. 9747839 | 
|| x () 76231164 | 
[a3] kK, 8 = 05 x @ 3.85285 | 
| >| Ys zs Vr - @ ; : } 145 9Olule5 
|@5| Ms FROM TABLES (@(4 | 2.807620 | 
46) | \ Mel = KS*! | 246234958 | 
[29] As/axfROM TABLES <5) | 3. 5255888 
[@@| SIN K/g = SIN (3 | 006719425 | 
[@| cos K/e = cos | 0699773991 | 
{@o| (6x (9) - | 25503722 
mE ( x (9) + © 2. 681.7607 
(2 9») + (16) X Ce) | 1. 1740237 
| 23 (9) - (© x ® | 0-8214561 | 
[| (x @& 8.9915637 | 
[@5! @ » el) 9- 4653623 _| 
KL) ©, @ be 1391249 | 
|e@)| «* 2 - At/a» -1. 3384536 
| 28) (2) x @4 +2 2. 8042165 
| (29) es: ly ee) | 0. 3566058 | 
[G9] 1 -9=1-@9 | 046433942 _| 
GD] €dxC) + GD XC —Ar/as | =0. 5106559 
32 rae = le yy | 0. 3815268 
| 33)| =1-@ __| _0.6184732 

SORORCHer @5)+ G3)xAr/xy?-l| 19. 3826585 
4 001745329 XB)xG4)+ (1)-{6)x7) 3. 4032603 
| 3 | __ 5729578 ra : 7 7+ 44355062 
Gp, 4 = G)xGQ Arfhx 5. 9758046 
| 38 | - = 62 +\/m2-1 8. 804237 











Mach Number and the physical characteristics of the 
wind tunnel. Similarly 6,* is used instead of 6,* (at 
point U) 

The value of 6,* may be estimated by methods such 
as that given in reference 6. By making the assump- 
tions that the pressure gradient is negligible, that the 
velocity-profile parameter V is constant, that the effec 
tive start of the linear boundary layer is at v4 /3 and the 
integration interval is from x4/3 to D, Eq. (7) of ref- 
erence 6 for a thermally insulated surface may be written 


in the form 


6p* = [0.051(3D — xa)I7)/(10 + Mr*)R“] (41a) 


where R is evaluated at the arithmetic mean of the 
static and wall temperatures in the test section and 
based upon the distance D — x,4/3, and // is obtained 
from Table VII of reference 6 for My and N = 2.2R' 

Considering the other assumptions made in any bound- 
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ary-layer calculations, this rather simplified treatment 
appears to give satisfactory results. 

Although the boundary layer is assumed to be zero at 
the throat and its rate of growth is also zero, its second 
derivative is 36,*/x,°so that it has a slight effect on the 
throat radius of curvature. 


Thus 

L/R’* = (/R*) + Cb" /2£) (42) 
so that 

R’* = R*/[1 + (86,*R*/x,*) | (43) 
If },; is the ordinate of the initial region with the 


boundary-layer correction included, then 


Yy;=y+ 6 


| (a — I)x4 tan 64 364* x \* 
=y+ + a 
2(n — 2) yA Ps (44) 
6,” ( % y x4 tan 64 (: 
2 Ns n(n — 2) \vy 


It may be noted in the above treatment that the con 
tour is corrected only for the boundary layer on the con 
toured walls. This type of correction will give more 
nearly parallel flow in planes normal to the contoured 
walls but will cause a small longitudinal pressure 
gradient to exist because of the boundary-layer growth 
on the flat walls. 


flat-wall boundary 


If is desired to correct the contour also 


for the layer, the 6* in the above 
equations should be the sum of the displacement thick- 
ness on the contoured wall plus that for one flat wall 
THIRD DERIVATIVES OF THE CONTOUR AT THE INFLEC 


TION POINT 


When the 
flexible-plate supersonic nozzle is fairly 


Mach Number range to be covered by a 
large, it be- 
comes impracticable to maintain the inflection point for 
each Mach Number at a given support point. It is 
that the 
initial and terminal curves be identical at the inflection 


therefore desirable third derivatives of the 


point. The third derivative of the analytic curve used 


for the initial region is relatively easy to determine. 


—Y,,'" = (364*/x,* [(m — 1) tan 64 x47] 


( 15 
boundary-layer correction in the terminal 


it does not affect the 
The curvature 


Since the 
region was assumed to be linear, 
third derivative of the terminal curve. 
is defined as 

C = d0,/dsy = (d*y/dx?) cos* 6, 16) 


The rate of change of curvature with s,, is 


dC/ds, = (d/dx)|(d?y/dx*) cos* 0,|(dx /dsw) | 
— (d*y dx*) cos! 6, — (47 
3(d?y/dx*) cos® 6. sin 0,.(d6, 7” 
also 





Howe 


and 


There 


or 


where 
deriv: 
sectio 
for th 
for tl 
which 
equal 
point. 
error. 
throat 
K 6, 
6; fou 
quiret 


In « 
precec 
herein 
of the 
perfec 
tion ¢ 
lor de 
of n, ¢ 
throat 
tive tl 
is the 
and tl 


Gen 
height 
bound 
proxin 


tment 


ETO at 
second 
on the 


(42) 


con 
> cOn- 
more 
yured 
ssure 
owth 
r also 
bove 
hick- 
wall 


by a 

be- 
it for 
It 1s 

the 
tion 
used 


le. 


}()) 





DESIGN OF TWO-DIMENSIONAL SUPERSONIC NOZZLES 


ac ds (d ds, \(d6, ds, ) = (d dt, (a8, ds, x 
(dt, dS, 
(d ‘dt,)(d6,, dt,) (dt, ds, 


d dt;)} [(@, — 0:) K](d0, dt); X (48 


(dt, ds, 
-(1SO0 7K )(d@, dt; (dt, ds)? 


(K in degrees 


However, 


(d0,/dt,;)? = 1/r,? (49 
and 
(dt,/dSy)a” = 1/(Ma? — 1) 50) 
Therefore 
(dy ‘dx*), = (dC/ds,) sec’ 04 (51) 
or 
i m7 sec! 0, 64(A7r/A*) |? : 
M4 7 ISOK(M,? — 1) |e A”) | eats 


(6, and A in degrees) 


where the subscript f is used to denote that this third 
derivative is for the terminal curve. For a given test 
section Mach Number, a given value of the exponent 1 
for the initial curve, and a given value of the ratio A 6,4 
for the terminal curve, there is only one value of 6, 
which will make the third derivative of the initial curve 
equal that of the terminal curve at the inflection 
point. This value of @, must be found by trial and 
error. If this value of 6, does not satisfy the length and 
throat radius requirements, then a different value of » or 
K 6, must be selected and the corresponding value of 
6, found to satisfy the continuous third derivative re- 


quirements. 


COMPUTING FORMS 


In order to simplify the computations described in the 
preceding section, computing forms are presented 
herein. Table 1 is for determining the terminal length 
of the nozzle (made nondimensional by dividing by the 
perfect fluid test-section half-height) for any combina- 


tion of 0; and A 64, (4v7 8; > S+ A 6;). Table 2 


is 
lor determining whether or not the chosen combination 
of n, @;, and K 6, will satisfy the conditions of length 
throat radius of curvature, and continuous third deriva- 
tive through the inflection point. Item 19 of this table 
is the difference between the physical half-height, //, 
and the perfect fluid half-height, 4: thus, 
bp* = II —h (53 
Generally, the test section Mach Number, half- 
height, and flexible-plate length are known, and the 
boundary-layer displacement thickness can be ap- 
proximated. Therefore, itis possible to enter Table 2 and 


Hu] 


solve for a combination of A’ 6,, 64, and m that produces 
a continuous curvature across the inflection point and 
then find the corresponding values of / h, R* h, Ax h, 
and Lh. If any of these values is not acceptable be 
cause of other considerations of the design, a new value 
of A 6, or n must be assumed and the problem resolved. 
It is to be noted, however, that the value of L 4 may be 
such that the nozzle terminates either upstream of, at, 
or downstream of the end of the flexible plate. Such a 
solution is valid and practical. All dimensions are in 


inches and all the angles are in degrees. 


APPLICATION 


As an illustration of the method, a hypothetical prob 
lem has been assumed and solved on Tables 1 and 2. 
For this problem a test section Mach Number of 3.0, a 
test section physical height of 40 in., and a flexible 
plate length of 200 in. were assumed. In addition, a 
boundary-layer thickness at the test section was as 
sumed to be 2.0 in. for the purpose of this illustration. 
Since the maximum value of A’ 6, is 4.0, this value was 
chosen together with a value of m = 3. These choices 
are arbitrary but convenient since a reduction in the 
value of A 6, will reduce the overall length, and an in- 
crease in the value of » will increase the overall length. 
If the final overall length is not satisfactory when the 
problem is solved for continuous third derivative at the 
inflection point, adjustments in these parameters are 
required. 

The initial value of 6, was chosen as the maximum 
value defined by the inequality stated in Eq. (12). 
Progressively smaller values were assumed until the 
curvature difference between the initial and final curves 
reached zero (Item 32, Table 2). Using this value of 
64, the terminal length was then computed on Table 1. 
Finally, the overall length and the radius of curvature 
at the throat were computed by completing Table 2. 


Mach NUMBER AND SECOND DERIVATIVE 
DISTRIBUTION ALONG CONTOUR 


The slope and second derivative of the initial curve 
can be obtained analytically. The values of Mach 
Number upstream of the inflection point can be ob 
tained approximately from the area ratio (before correc 
tion for boundary layer). Instead of using the actual 
ordinate at these points to obtain the area, the distance 
along an are normal to the profile at each point is sug- 
gested. This procedure makes the Mach Number at 
the inflection point consistent with the value calculated 
at this point for the terminal curve. 

For the terminal curve, the Mach Number can be ob- 
tained from the Prandtl-Meyer expansion angle 


Ve = va + On — I | 
= pr — Oy — (K /4) — (04 — 8, + (54) 
V (04 — 0,)2K 


from the inflection point A to the characteristic point P 


and 
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TABLE 2 
Nozzle Lengths, Throat Radius of Curvature, and Third Derivative 
Mr = 3.0, vr = 49.7573, Ar/A* = 4.23457, H = 20,1 = 163.7129, 3D = 600, VW Mr? — 1 = 2.828427 
1 | & 16. 58577 12. 9. 7.0000 7.7075 77057 | 
2 | kK‘, lhe he he i he bh» | 
 uaensor eee einer - ———_ 
3 | K= Ox® 66+ 34308 48.0 36. 28. 30. 83 30. 8228 
4 |K/4 = @4 16. 58577 Cs a a 7.075 7.7057 
5 |y,. = = -@)-@ 16. 58577 25-7573 | 31-7573 | 35-7573 | 3a 3423 Be 3459 | 
6 |Ma FROM TABLES@Q(G)| 1.658399 1.977493 24-2009 || 26 3594K6 =| 2 30239 2. 302518 
7 | AsA*FROM TABLES@6)| 1.299611 1. 65630 200658 _ i 2.31521 | 219791 | 2.198143 
8 [sn O,= SING) 0. 2852845 _0.2079117 0. 156434 001218693 | 02134116 | 0.134087 
9 | TANG, = TAN © ___|  o.297emi9 | 0.2125566 | 0.158364 | 0.122785 | 0135339 | _0..13530661 | 
10 | 0.01745229 x @ 0. 2894763 | 0. 2094395 | _0. 1570796 0. 122173 O.134521 | 0.13448982 | \ 
iu 2b oO 1. 280792 1. 6442177 1. 998332 | _ 2+ 309455 2.191292 2.1915227_ | hie 
te : a oe 3- eee SS eae SS super 
2 /Ar | 
13 a diem O- 3542272 O- 3542272 O- 3542272 | 0. 3542272 03542272 | 0» 35422722 sists 
14 |G) - 0. 280792 Oe 64442177 0.998332 | 1.309455 1.191292 | 1.195227 minit 
- a f yresst 
15 | Xa/h = G@x@/@ 0. 339469 1.07359% 2.232778 | 3.777697 3.118008 | _3-1193581__ me 
16 | 4% = BP | 5° 9758046 | carrie 
=~ = T 
17| %, =©+© | | 9..0951627 ized 
iA — 
18 | h = L @ 18.0 18.0 ‘ss 18.0 | 18.0 18.0 | 18.0 boun 
19 a= H -(8) 2.0 2.0 | 20 | 2.0 2.0 | 20 = 
nities a ee — : as ni 
20 | Xa =@x 6.011044 | 19.324710 |  40.190004 | 67998546 | 5622414 | 561484458 volvi 
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Vo = vr — 0, (55 and from point 27 + 1 to point 27 + k 
from the characteristic point to the end of the nozzle. 6/04 = (4r + 2k + 1 — 2m)?/16rk (59 s 
The value of 6,, can be readily obtained at the points ob- ; ; 
9 shane ite oa ie aes where k, m, and r are all integers. Although the inflec- “- 
tained by automatic computing equipment. For this : v,. 
oe ie a ee as ; tion point A is usually given the number zero, according s 
computation, the are AO (Fig. 2) is divided into 7 equal é i ; : 
ee ie : 5 to the above notation, it would correspond to a point 3 
parts. The angle A 2 contains k of these 7 parts so that nae : : cao as 
number of |». Likewise, the characteristic point | 
kr = 1 — (O:p/64) = K/204 (56) would correspond to a point number of k + ! 2, the Y 
= ors. ; <a int Q (Fig. 2) would correspond to a point number of oo 
The characteristic network is established by initiating - iy be nae nig Aetecehies -nitens = é 
rte oes ; : 2r + 1/2, and the point U (Fig. 2) would correspond toa 
the characteristics on the are AO at the midpoint of : ' 
z ae ‘ point number of 27 + k + 1)». ¥ 
each of the y parts. The points on the contour are num- a : : , : ’ 
: a ; ‘a : : The slope of the terminal curve is tan @,, before the 
bered successively from | to (27 + k), point 1 being the tees ; ' ‘ 
a : . j correction for boundary layer is applied. After correct- Rec 
first point following the inflection point. Let m be the . ole ar es 
; : . : ing, the slope of the contour is tan 6, + 364* 2x4. Using .?T 
point number; then, from point | to point , ; AE 
the values of tan 6,, as computed above together with the | A Ep 
6,,/04 = 1 — [(2m — 1)?/16rk] (57) values of x as obtained from the automatic computing Engin 
, 2 : is ar orm eA SIO > ob- operat 
trom point & + 1 to point 2y, equipment, values of the second derivative can be ob +c] 
tained by the method of divided differences. t Pr 
** E 


6./04 = 1 — [(2m — 1 — k)/4r] (58) (Continued on page 732) 
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Aerodynamic Flutter Derivatives for a 
Flexible Wing with Supersonic and 
Subsonic Edges’ 


SAMUEL PINES,* JOHN DUGUNDJI,? anp JOSEPH NEURINGER** 


Republic Aviation Corporation 


SUMMARY 


A box method is developed for obtaining the generalized air 
forces on an oscillating flexible wing in supersonic flow with both 
supersonic and subsonic edges. Essentially, the method con 
sists of representing the wing by a grid of square boxes and deter- 
mining the influence of one box on another These aerodynamic 
pressure influence coefficients when tabulated, permit the flutter 
analysis of an arbitrary wing with arbitrary normal modes to be 
carried out in a routine way. The coefficients satisfy the linear 
ized unsteady supersonic flow equations and the downwash 
boundary conditions, and, in addition, are formulated in a man- 
ner independent of the modal shapes of the structure 

The box method is applied to some simplified examples in- 
volving rigid body modes, and agreement with other methods 
is seen to be reasonably good 

The box procedure appears to offer a simple routine manner 
of analyzing flexible wings for supersonic flutter analyses which 
is well adapted to programming on computing machinery. How- 
ever, the square box method as developed here for subsonic edges 
is inapplicable to Mach Numbers below ./ = 1.414 without 
further modification. Some suggested modifications for extend 


ing below the J = 1.414 range are discussed in the body of the 


paper 
SYMBOLS 
speed of sound 
= vertical deflection of airfoil positiv e down 
k = reduced frequency (ws/ 
VW = free-stream Mach Number 
p = pressure 
R = hyperbolic distance (V(x — 4)? - 3*(y — n) ) 
R il = pressure influence coefficient 
= length of square box side 
= time 
= velocity 
V = area included in forward Mach cone 
= downwash (positive down 
a F = coordinates 
ee a = nondimensional coordinates 
a = angle of attack 
3 = the definition (1/ M2 — 1) 
= the definition (.17?/1/* — 1) 
9 = angle of roll 
5.7 = coordinates 
9 = nondimensional coordinates 
p = fluid density 
¢ = velocity potential 
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w = frequency of oscillation 


: hin w M? 
w = the definition ) 
U MW l 


= column matrix 
= square matrix 
row matrix 
Subscripts 


= center of box 


D = diaphragm 

j = any box 

,n = box orientation subscripts 
HW’ = wing 


INTRODUCTION 


gage PROCEDURES have been devised for obtain- 
ing the supersonic three-dimensional flutter de- 
rivatives for wings oscillating in elastic normal modes. 
References 1-5 are good illustrations of some of the 
various approaches to this problem. In general, these 
approaches have limitations either as to type of plan 
form—they apply to delta wings; as to type of deflec- 
tion shape—they apply to specified linear, parabolic, 
or cubic mode shapes; or as to frequency parameter 
they apply only to low values of reduced frequency pa- 
rameter, k. 

In order to avoid these limitations and to attempt to 
obtain a more general method for arbitrary wings with 
arbitrary deflection shapes, a numerical integration 
procedure was developed at Republic Aviation Cor- 
poration. The method consists of representing the 
wing by a grid of square boxes (see Fig. 1) and deter- 
mining the pressure at any one box caused by the 
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Fic. 1. Box representation of a typical delta wing 
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motion of all boxes in its forward Mach cone of de- 
pendence. This requires the evaluation of a set of 
pressure influence coefficients which will be functions 
of box location, Mach Number .\/, and reduced fre- 
quency parameter k. The box procedure has the ad- 
vantage of requiring an accurate representation of the 
integral kernel over a much smaller portion of the 
wing surface, thus permitting expansion procedures to 
hold over a large reduced frequency range. In addition, 
it becomes possible to separate out the elastic defor- 
mation motions from the integrals, thus permitting 
the evaluation of pressure influence coefficients which 
do not depend on the wing mode shape. Wings of 
arbitrary plan form can be treated, since the box 
method, as will be shown, is readily adaptable to sub- 
sonic as well as supersonic edges. 

A box procedure was originally developed by V. 
Zadikov and J. Neuringer as a means of dealing with 
static aeroelastic effects on supersonic wings.’ The 
procedure permitted the determination of the effects 
of chordwise and spanwise bending on loading pat- 
terns. It was subsequently adapted to unsteady mo- 
tion using velocity potential influence coefficients.’ 
Further refinement of the unsteady motion problem, 
including change to pressure influence coefficients and 
adaptation to subsonic edges, led to the present form 
of solution.* * The material described in this paper is 
essentially a presentation of the work reported in refer- 
ences 8 and 9. 

DEFLECTIONS ON THE AIRFOIL 

The wing is divided into a lattice of square boxes 
which completely cover the surface as shown in Fig. 1. 
Each wing square is assumed to be oscillating har- 
monically as a flat plate about its center (x,, y,) with a 
vertical deflection given by 
h(x, y, 0) = et Jhi(xe, Vo) + 

(x — x )aj(Xe Vo) + (Y — Ve)O(Xe Ve (1) 

The instantaneous vertical velocity at each square 


box is given by 


oh _Oh 
w= + [ = 
Or Ox 
eV iwh; + iw(x — x,)a; + iwl(y — y.)0; + Ua (2) 


The contributions to the lift due to terms with 
(x — x,) and (y — y,) tend to zero with smaller box 
size, so that these terms may be neglected when a wing 
plan form is represented by many small boxes. For 
this analysis then, the final form for the vertical ve- 
locity at each box becomes 


tf. , ‘ 
w, =e” \twh)(X Vo) + Uaj(x,, y,)} (3) 


Thus, each wing square has two degrees of freedom, 
viz., translation and local angle of attack. 


AERODYNAMIC EQUATIONS 


The equations of motion for unsteady compressible 
potential flow, when linearized, can be expressed by 
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the single partial differential equation 
2M Oe | Oy 


O-¢ O-¢ O-¢ oe 
Ox Ot c* Of 


(1 — JT) — 
Ox* Oy* Os° ( 


In the above equation, the dependent variable 
represents the perturbation velocity potential, and 
AJ represents the free stream Mach Number of the air 
flow. <A solution of Eq. (4) which, for the case of an 
airfoil with supersonic edges, also satisfies the harmoni- 
cally oscillating downwash boundary conditions on the 


z = Oplane is given by" 


—e"(* "w(é, nje ‘**~S’cos (@/M)R 
I dé dy (5 


R= V(x — §)? — Biv — 9 


where w(é, 7) represents the downwash velocity on the 
airfoil and where the integration region |” extends over 
the forward Mach cone region of any receiving point 
(x, y). See Fig. 1. Since the downwash velocity w 
is a varying function over |” but assumed constant in 
each box [see Eq. (3)], Eq. (5) may be written as 


— eit i 
g(x, y, t) = >» W; (Xj, Vj) X 
7 
l = * cos (@/.)R 
J Vie R 
where |’; is the area of box 7 included inside of the for- 


ward Mach cone emanating from point (x, y). 
The pressure at (x, y) is given by the expression 


¢ O¢ 
p(x, y, t) = 2» ( + U ) (7 
; Of Ov 


Upon substitution of the velocity potential as given 
in Eq. (6) into Eq. (7), the pressure at point (., y) be 


didn (6 


comes 


~ Dae .. Uo 
p(x, y) = ee > w|i + | ‘ 
" w Or] 
PP gt 0 on (2/18) B 
dé dn (S 
J Vin R 


At this stage, it becomes convenient to nondimension 
alize the coordinates with respect to the length of the 


square box side, s. 


Let 
i. é = x's. t's 
y, n = V's; n’S | (9 
ws/U =k | 
M?/M?—1= y 


The final expression for the pressure at any port 


, 


(x’, y’) then becomes 


} 


p(x’, v’) = 2pws e' DY (R; + iw (10 
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where the pressure influence coefficient, R; + 7/; is de- 
fined by 

-i.,19 

ae 1+ : 

R k Ox 


us 


x 

a ee , 
,—tky ) aos (by/M) 

| é = (ky/M) R dt’ dn’ (11) 

Vie R 


The pressure influence coefficient as given in Eq. (11) 
is representative of the oscillating pressure at point 
v’, y’) caused by an oscillating unit normal velocity 
of the box V’ 
seen to be functions of Mach Number /, reduced 
frequency k, and the differences of the coordinates 
between the receiving point (x’, y’) and the grid box 
!.. This dependence on differences in coordinates 


The pressure influence coefficients are 


can be conveniently altered to some orientation identi- 
fication scheme. Such an orientation scheme for a 
square box grid system is indicated in Fig. 2. For ap- 
plication to a wing plan form at a given ./ and &, it is 
necessary to have available a table of the numerical 
values of R + c/ for different orientations. This single 
table can then be used for determining the total pres- 
sures at all grid boxes. 

The application of the pressure influence coefficients 
to the determination of the flutter derivatives follows 
readily from Eq. (10). Rewriting this pressure equa- 
tion in matrix form, 


1p! = 2pwse’™ (R+ 11) X jw} (12) 


and introducing into the above the expression for down- 


wash over the airfoil, Eq. (3), 


h l 
rf = ‘6 J { 73 y¢ ' ( 2) 
»>wy = w ar rw b ) ay (1: 
vields 
er fh — a 
2 pw°s*e"™ (—/+7R) >+ (R+i1),;ay>> (14) 
| ish ok \\ 


The work done by these aerodynamic forces is given 
by 


aie | h(x, vy) p (x, y) dx dy (15) 
e wine 


The assumption is now made that the pressures are 
constant over each box. Eq. (15) can then be combined 
with Eq. (14) to give the following matrix form for the 
work done by the pth mode against the aerodynamic 


forces of the gth mode: 
W oe = 2 pwsre'™ x 


h, ha . 1[h, | 
"r+ init y! *1R + iD) fa} )) (16) 


Eq. (16) represents the final form for the flutter de- 
‘ivatives. It is seen to be in convenient matrix form 


and to involve matrix multiplications of the mode 
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Box Tree = (hn) 
THUS BOx H HAS L-3,N-4 
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Fic. 2. Notation scheme for pressure influence coefficients 
shapes of the wing with the pressure influence co 


efficients. 


PRESSURE INFLUENCE COEFFICIENTS 


The pressure influence coeflicients, R + i/ defined 
by Eq. (11), have been evaluated by making use of 
some approximations in treating the exponential and 
cosine factors inside the integral. For boxes in the 
immediate vicinity of the receiving point (values of 
n < 1), an expansion of the exponential and cosine 
factors up to terms in k? has been used. For boxes 
further removed from the receiving point (values of 
n = 2), the mean values of the exponential and cosine 
factors were used. It should be remarked that in ob 
taining the portion of the pressure coefficient arising 
from the O Ox’ term, the differentiation with respect 
to x’ was performed before making approximations on 
the exponential and cosine factors inside the integral. 

In terms of the orientation identification scheme indi 
cated in Fig. 2, the pressure influence coefficient R,, + 


i/,, can be written as 


R ’ . ll 10 / 4 l 
w tiln= lit pal, J 
; k> 
e ** cos pe VX? - a7) . | ’ 
-- —— dX d¥Y (ij 
Vv A? — BY: 


The above formula applies specifically for a square 
which lies completely within the forward Mach cone of 
the receiving point. This formula can also be made to 
apply for any other type of box whose edge is cut by 
the forward Mach cone provided the definitions as 
stated in Tables 1 and 2 are observed. For the case 
of n = 0, the lower limit » — (1 2) should be replaced 
by zero. 

The algebraic evaluation of Eq. (17) for the pressure 
coefficients is given in Tables 1 and 2. These results 
were obtained using the approximations to the ex 
ponential and cosine factors described previously.* 
A more exact expression, which does not make any 


* The authors wish to thank Messrs. Holt Ashley and G 
Zartarian of the Massachusetts Institute of Technology for 
pointing out a numerical error in one of the coefficients of the 
original expressions appearing in reference 8 
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TABLE 1 TABLE 3 : appr 
Formula for Approximate Pressure Influence Coefficients Table of Pressure Influence Coefficients for 17 = 1.1414, k = 0.02 Wath 
(m = 0, 1) : 
This 
M =L4y4 h=.02 
(be. =| -/B(f- B(ht. : e invol 
ri daft BP <i A rtm le} = 
. * i. r “IB{ De of - oO . 
+h fer sieos" AED. n-sheos HOS nf -5feos RU) tf " R be . 
| for 
ie 1 B(L+,5)) , aN ry. 2. A2/ Pp - 5) slwa’ 
nf -s}eos n-5 | tay l(t 5) \/(n+.5)*-B°(L-.5) eo ti 
= hence 
~(£+.5) nes) B°(b45)* -n(£-.5) in- SFOS? lle sPOa} As 
influe 
Sealy -1N+.5 _ HIN$.5 no =I -.5 two-( 
I] ak 5)Cosh B(E-5) (b+.5)cosh Brass) n(g S)cosh B(L-.5) be cé 
-~+n-. -iB(L-.5) of ' 
+n(b+5)cosh G Say} al ine sdes nes roe 
( 
~(n+5 Jeog!Bbeg) -nn- Sco Mesos EI} two- 
THE FOLLOWING DEFIN/TIONS HOLD IN THE ABOVE EQUATIONS R 
cos ‘x=0 FoR x >+/ 
cos’x=0 toe FoR O£ XK£/ 
COS'%=T-coS tx) FOR ¥< 0 
vx = 0 For x<0O 
cosh x=Qto +00 FoR X>+/ 
cosh'x= 0 FoR 0£%£+/ 
coshx=coshex) For x<0 In tl 
ordet 
® wh 
In 
flutte 
effici 
TABLE 2 cues 
Formula for Approximate Pressure Influence Coefficients -8. 8F0| & 6¥20 lrequ 
(nm = 2) 
, al 
R--4 sin(n &¥) cos (A Vn @L) X {B} TI 
/ RE 7 sonic 
- = costin+.F 0 Vin+5)2- ig? x}A +} a : 2 
Ret 5([n+5].Rx)cos("y BC) x} An M=I4I4 ¢-1 -.02 h,* .075 KX, 0 vie 
he Y, be tt 
+ <4 cos([n-.5 05 Gein P)x An- ' 
ROT s(( Lh ¥)< M 6 | from 
if (\ subst 
I- - £cos(nht)cos(A vnt=Bee) {8} yy phra; 
ihe —_ wing 
+ pa sin({n +.5]h¥)cos (A Kins -F2*) x { Ane n/h a2 | 2 4\n | dai 
a " Pe cK 7 I sub 
, 32 hy | 3 | 3 | 4 32 | — 
~ pepsin]. a) 0s fet N-SIRRE) fo} . 3. diap 
w | 3/R6| 5 | 5 | 6 439 | 94 nei 
WHERE — aa Lt d 7” 4— - , 
, = wo} ar ha| 817) 7 8194 3 efficie 
Pe nt. Ny N+.§ =! Tie, —— = "e 
{B} = -(f-5)cosh bas +(2+.5)cosh aL +(2-4)cosh 6-5) siwiwtelslelelslepe this | 
N usua 
: AT 
~(04.5)cosh =F 4 M48 _ne y2|\4l4o hw) | 4) | 3 | |e | oN —— 
(Bes)cosh aires) B Ant} "at | An- : oo 
. ¥ ws | 44 | 43/820] 19 | IB | 17 | 17 | 18 | IT | 20 — 
f Ane} = cos’ S22). cos" (L+.5) | sage tribu 
nt.5 nt. 46 y | 24 | 23 | 22 | 20 | a | 22 | 23 | 24 — 
{a = cos BCL-.5) _ cos! B(L+5) =a 30 l 29.| 29 | 27 | 26 26 27 [as 29 Als 
n- N-.5 ne-.5 ZL. Bassvaa RRR CRNRS CRNUNY 2/208 DAE M) FHT aE z 
y | illust 
| scher 
THE SAME DEFINITIONS HOLD AS FOR THE — os 


n=0,/ FORMULAS Fic. 3. Subsonic leading edge wing layout sti 
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:pproximation on the integrand, was worked out by 
Watkins of the NACA" and is given in reference 8. 
[his expression is somewhat more complicated as it 
functions and are 


involves infinite series of Bessel 


sines. It is felt that over the range of frequencies to 
be used for most flutter problems, the approximate 
formulation is adequate. In any event, accuracy can 
always be improved by going to a smaller box size and 
hence lower reduced frequency R. 

As a check on the values computed, the pressure 
influence coeflicients can be combined in groups to form 
two-dimensional triangular pieces. These groups can 
be compared with the exact two-dimensional values 
of pressure for constant unit downwash.!” For a two- 
'/5)s, the exact 


dimensional strip of length x = (m + 


two-dimensional pressure is 


—] (=) | MW- 1, 
Je cos w — w (fo), 
kp M M? 


—] > aa MP |... 
I = ba — Ji VV sin w + We w(fo)r 


= 


(1S) 


In the above, 
order, and (fy) and (fo), are functions of both 17 and 


J) represents Bessel function of zero 


® whose numerical tabulation is given in reference 12. 

In any practical application of the box method to 
flutter, it is assumed that tables of these pressure co- 
A typical set of pressure co- 
1.414 and reduced 


efficients are available.* 
eficients for Mach Number J/ = 
frequency k = 0.02 is given in Table 3. 


TREATMENT OF SUBSONIC LEADING EDGES 


The preceding work applies only to wings with super- 
sonic edges since the basic solution for ¢ given in Eq. (4) 
can satisfy the boundary conditions only if the top and 
bottom surface of the wing are completely isolated 
from one another. To treat the case of a wing with a 
subsonic leading edge, the Evvard'* concept of a dia- 
phragm of unknown downwash is added to the given 
wing. This diaphragm, which covers the region be- 
tween the foremost Mach line of the wing and the 
subsonic leading edge, results in a combined wing plus 
diaphragm plan form area which has completely super- 
sonic edges (see Fig. 3). The standard pressure co- 
efficients for supersonic edges can then be applied to 
this combined wing plus diaphragm plan form in the 
usualmanner. The downwash of the diaphragm, which 
is determined by applying the condition of zero pres- 
sure at a diaphragm point, creates an additional con- 
tribution to the pressure at any point on the wing sur- 
face. 

Algebraically, the treatment of subsonic edges can be 
illustrated in the following manner. The basic ordering 
scheme for the pressure at any point on the wing- 


* Wright Air Development Center is presently computing 


some of these pressure coefficients 
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diaphragm regions is given in matrix notation by 


jP\ _ | pwn pu D jw 19) 
lof Pow ‘Pop Lewy | , 
where 

P = total pressure on wing box 

w = known downwash on wing box 

Wp = unknown downwash on diaphragm box 

pwn pressure coefficients on wing due to boxes 
on wing 

Pwo = pressure coefficients on wing due to boxes 
on diaphragm 

Pow = pressure coefficients on diaphragm due to 
boxes on wing 

Pop = pressure coefficients on diaphragm due to 


boxes on diaphragm 
5 


In the above, P, w, and wp are column matrices while 
Pbww, Pwo, Pow, and ppp are rectangular matrices. 

Solving for wp from the lower set of the matrix equa 
tions and inserting into the upper matrix equations 
results in the following expression for the total pressure 


P at any point on the wing. 


P} = (pwn Pwo Ppp" pow) X By 20) 
Thus, the effect of a subsonic edge is to introduce a 
correction term, PwpPpp Pow, subtracted 
from the basic wing pressure influence coefficients. It 
should be borne in mind that the elements of the p ma- 
Also, 


+14, 


which is 


trices are complex numbers. Ppp will always be 


a triangular matrix for JJ > 1. and hence its in- 


version is readily accomplished. 
APPLICATIONS TO SPECIFIC WINGS 


Subsonic Edge Delta Wing 
applied to a subsonic edge 
a rigid translation mode for 


box method was 


The 


delta wing oscillating in 


M = 1.414 and k = 0.02 (note that & is based on box 
size, § The leading edge angle of sweep for the delta 
was 63.4 The box representation of the wing to- 


9 


gether with the required diaphragm is given in Fig. 3. 
The basic ordering matrix for the pressure at any box 
The final lift and moment (about 
The 
box method values are compared with the theoretical 


values as obtained by Watkins’ and consequently, 


is given in Fig. 4. 
the wing apex) distributions are plotted in Fig. 5. 


his notation of lift and moment is used—namely, 
Lift = P = —4pbl kpc” X 
hy 
(L,; + tL2)— + (Ls + tL4) a 
b 
Moment = VM, = —4pb?U°k,*e" X 
hy 
(VU, + 1\Me) . + (lJ, + iV, a 2] 
D 


where in the above equations 6 is the root semichord 


of the wing and kz is equal to wh U 
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It is seen that the agreement between the two meth- 
ods is reasonably good for this case. For interest, the 
iift and moment as obtained from the uncorrected pres- 
sures is also plotted. This leads to a sizeable difference 
in lift distribution, particularly for the ZL; (imaginary 


part of the pressure coefficients). 


Sonic Edge Delta Wing 

The box method was also applied to a sonic edge 
limiting case of supersonic edge) delta wing oscillating 
in both rigid translation and rigid pitch modes about 
the apex. The Mach Number and reduced frequency 
were again taken as JJ = 1.414 and k = 0.02, and the 
leading edge angle of sweep was 45 The representa- 
tion of this wing was made up from the previous sub- 
sonic case by treating the former diaphragm boxes 
as wing surface boxes. A seven box root chord was 
used resulting in a total of 28 boxes over one-half of 
the wing. The final lift and moment (about the wing 
apex) distributions are plotted in Figs. 6 and 7, and are 
compared with the theoretical values from reference 
|. The agreement for these cases is seen to be rather 
good. 

The results presented here are described more fully 
in reference 9. Also presented in there is the ap- 
plication of the method to a rectangular wing which 
again yielded good agreement with theory. The pres- 
sures at each box as obtained by the box method are 
also shown in reference 9. From observation of these 
pressure distributions, it appears that the jagged 
leading edge causes oscillations in pressure which 
propagate themselves down along the Mach lines. 
Hence, the wing must be covered with a sufficient num- 
ber of boxes to average out these pressure variations 


caused by this jagged edge. 


DISCUSSION 


The application of the box method to determining the 
flutter derivatives for flutter analyses follows readily. 
The wing plan form is approximated by a grid of square 
boxes, the basic ordering matrix is determined for the 
lowest Mach Number, J = 1.414, and the deflection 
and slope at each box for each mode is found. Assum- 
ing the existence of a table of pressure influence co- 
eflicients, the remaining computation, although involv- 
ing many multiplications, is very routine and is well 
adapted for programming in a digital computer. All 
reduced frequencies k, and all Mach Numbers J/ 
greater than 1.414, can be done without altering either 
the box representation of the wing, the deflections and 
slopes, or the basic ordering matrix. Hence, a mini- 
mum of setting up work is passed on to the flutter 
analyst. Since subsonic edges can also be treated in 
the manner described in a previous section, Treatment 
f Subsonic Leading Edges, and programmed on a 
digital computer, the method is applicable to all plan 
forms for M greater than 1.414. 


{ 


The accuracy of the box method, within the frame- 
work of linearized aerodynamic theory, increases with 


the number of boxes used. <A large number of boxes 
seems necessary to average out the pressure variations 
caused by jagged leading edges. Based on the examples 
here presented, a minimum of about 30 boxes on one 
half the wing is required. For flutter analyses where 
tip effects assume greater importance, somewhat more 
boxes on the wing would be necessary. It should be 
mentioned that probably simple and routine corrections 
can be developed to account for nonsquare boxes at 
the leading edge of the wing, thereby refining the 
treatment of wing leading edges and reducing the 
jagged edge pressure variations. 

The use of square boxes leads to some difliculty when 
applied to Mach Numbers lower than 1.414. Below 
1.414, the forward Mach cone emanating from the 
center of a box cuts across the two sides as well as the 
front of the square box. Therefore neighboring boxes 
at the same chordwise location will influence one 
another (there is in addition to the 00 pressure con- 
tribution also 10, 20, 30, etc. contributions). This 
presents difficulty in applying the standard square box 
method technique to a box where a side boundary (either 
wing or diaphragm) ends, since beyond the side bound- 
ary, the pressure will not be zero. There is particular 
difficulty in the case of a subsonic leading edge, for the 
diaphragm must extend itself sideways beyond its 
correct boundaries towards infinity in order to insure 
that the pressure is zero everywhere off of the wing. 
The standard square box method technique therefore 
requires some modification to deal with the range 
below \/J = 1.414. 
ify the shape of edge boxes by clipping off a corner 
Another approach would be 


One approach would be to mod- 


or something similar. 
to make use of rectangular boxes below \/J = 1.414 
and thus prevent side neighboring boxes from influ- 
encing one another. This would bring down the limit- 
ing Mach Number of 1.414 to some lower Mach Number 
based on the rectangular box’s dimensions. Still an- 
other approach is to utilize diamond shaped boxes 
whose edges lie along characteristic Mach lines. As 
yet, these approaches have not been fuily explored. 


CONCLUSION 


A box method has been presented for the determi- 
nation of the generalized air forces on oscillating super- 
sonic wings with both supersonic and subsonic edges. 

Application of this general method to some simplified 
examples shows reasonably good agreement with other 
more limited theoretical methods. 

The box procedure appears to offer a simple and 
routine manner of analyzing flexible wings for super- 
sonic flutter analyses which is well adapted to pro- 
gramming on digital computing machinery. Accuracy 
within the framework of linearized theory can be in- 
creased by using more boxes. 

Some modification to the present box procedure using 
square boxes is necessary to extend down below the 
M = 1.414 range. 
gested but they have not as yet been fully explored. 


Several proposals have been sug- 
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An Analytic Method for the Design of 


Two-Dimensional Asymmetric Nozzles’ 


HANS PETER LIEPMAN?t 
Unwersity of Michigan 


SUMMARY 


An analytic method is presented for the design of two-dimen 
sional asymmetric curved nozzles. Such nozzles may be used 
to continuously vary the Mach Number by relative translation 
of one contour with respect to the other 

The continuity equation and the condition of irrotationality 
for steady compressible flow are given in curvilinear coordinates 
with stream function and velocity potential as independent vari 
ibles. Using Bernoulli's equation for isentropic gas and a 
suitably chosen analytic velocity distribution along the initial 
expansion contour, which is specified in terms of its curvature, a 
set of partial differential equations is obtained. These equations 
ire solved by expanding the dependent variables—velocity gq, 
flow direction @, and curvilinear coordinates (x, y) of a stream 
line—in powers of the stream function. This solution determines 
the nozzle contours for any constant value of the stream function 
as well as all flow characteristics in the subsonic region and the 
supersonic expansion region of the nozzle. From the general 
solution and the boundary conditions at the sonic line it follows 
that the curvature of the contour must be zero at the throat 

After the patching Mach line is calculated numerically the re 
compression contour can be determined to complete the nozzle 
calculations. These equations show that the Mach Number 
gradient must go to zero at the design Mach wave. This require 
ment implies continuous wall curvature at the corners of the 
recompression region 

The variation of Mach Number by relative translation of 
fixed contours of the asymmetric type is then discussed. An 

ration method is derived for this purpose which should con 
verge very rapidly as compared to the iteration required by the 


characteristics method of design 


NOTATION 
A = channel area (one-dimensional theory 
C; = aerodynamic force coefticient 
= speed of sound 
H = stagnation pressure 
h = area ratio parameter 
J = Jacobian 


= length of left running characteristic in recom 


pression region 


WV = Mach Number 

p = static pressure 

q = resultant velocit, 

diy 92, 9 = coefficients of g-series 

R = radius of curvature of reference contour 
S = are length along streamline 


= test section height 
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u = velocity component in curvilinear x-direction 
= velocity component in curvilinear y-direction 
Y = Cartesian abscissa 
y = Cartesian ordinate 
r = curvilinear abscissa (along reference contour 
“i, Xe, 2 = coefficients of x-series 
y = curvilinear ordinate (normal to reference contour 
v1, ¥ ¥s = coefficients of y-series 
a = Mach angle 
7 = ratio of spectfied heats 
6 = ho-perturbation function 
é = angle between neighboring radii of curvature 
n = independent variable, proportional to stream 


function 

nN = design value of 7 on upper contour 

7] = flow direction measured from curvilinear abscissa 
(on reference contour 

6;, Oo, Os = coefficients of @-series 

kK = curvature (= 1/R 

= independent variable, proportional to potential 


function 


p = density of gas 

¥ = stream function 

Y = potential function 

w = slope of reference contour measured from Car 


tesian abscissa 


Subscript 


0 = conditions along reference contour 
D = test section design value 
l, u = condition on lower, upper contour, respectively 


Superscript 


= critical (sonic) state 
INTRODUCTION 
The Asymmetric Curved Nozzle—Principle and Advantages 
p | (aig TYPE OF ASYMMETRIC curved nozzle under con 
sideration is shown schematically in Fig. 1. Re- 
cent experimental work by the NACA! * has indicated 
that such a nozzle can be used to vary continuously 


the test section Mach Number by the relative transla- 


tion of one contour with respect to the other. The 








Fic. 1. Sketch of asymmetric curved nozzle 
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Fic. 2. Geometry in curvilinear coordinates 


University of Michigan is presently carrying out a 
sponsored research project* directed toward the ex- 
tension of the Mach Number range of such a nozzle 
and toward the experimental determination of the 
nozzle contours resulting in most uniform test section 
flow over a nominal range of Mach Numbers from 1.4 
to 4.0.° 

The mechanical simplicity of translating one contour 
block with respect to the other to achieve a continuous 
variation in Mach Number is the most outstanding 
feature of this nozzle. Furthermore, potentialities exist 
with such a nozzle to simulate accelerated flight condi- 
tions by time-programming the nozzle translation. 


Design Methods for Asymmetric Nozzle Contours 


At the present time only the graphical method of 
characteristics is available for the design of the contours 
of such an asymmetrically curved nozzle. This method 
generally requires that the shape of the sonic line be 
assumed initially before the characteristic diagram can 
be started. 
the only guidance available is the one-dimensional 


For the design of the subsonic contour 


area theory which can be interpreted rather freely, 
especially in a curved channel. 

Experience with symmetrical nozzles has indicated 
that an analytic method can be of considerable value 
in certain nozzle design problems. Similarly, an an- 
alytic method of design for asymmetric nozzles might 
supplement the method of characteristics. An analytic 
design method was therefore developed from the most 


* The sponsorship of this project by the Aeronautical Research 
Laboratory, Air Research and Development Command, U.S 
Air Force, is gratefully acknowledged. Initial work of the proj- 
ect showed the lack of an analytic method of design which di- 
rectly motivated the development of the theory presented in this 
paper as an academic assignment supplementary to the spon- 


sored project. 
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suitable nozzle design procedure for symmetric nozzles 4 
This method does not require any assumption regarding 
the shape of the sonic line and extends into the subsonj 
region. 

Neither the method of characteristics nor the ney 
analytic design method achieves the variable Mach 
Number feature without a process of iteration. The 
iteration for the analytic design method, however, 
appears to be considerably shorter than that for the 
characteristics method. Of course the high degree of 
numerical accuracy of the analytic method need not be 
lost in this iteration process. 


PERFECT FLUID THEORY OF ASYMMETRIC NozzLe 
CONTOURS FOR A GIVEN Exit MACH NUMBER 


Fundamental Equations in Curvilinear Coordinates 


In view of the curved asymmetric contours, the use 
of curvilinear coordinates suggests itself to describe 
the flow field in regions I and II of Fig. 1. These 
coordinates may conveniently be chosen along and nor- 
mal to a reference curve Cj, usually the lower contour. 
Let (x, y) denote such curvilinear coordinates and let 
(X, VY) stand for the Cartesian base. The relations 
between (x, y) and (X, )), as well as the curvature « of 
the reference curve, can be derived from Fig. 2 giving 


X = Xo + [ cos w dx — y sin e| 
/ 0 


(1) 
Y= ¥Yo+ f sin w dx + y cos w | 
0 
k(x) = 1/R = dw/dx = (d?Y/dX?)/{1 + | 
(dV /dxX)?)}*? (2) 


w =Sx(x)dx = tan-! (dV /dX) = tan! F’(X) 


For steady irrotational flow of a compressible isen- 
tropic gas, the equation of continuity and the condition 
of irrotationality are satisfied by defining a curvilinear 
stream function and velocity potential ¢ as follows’ 


Oy /Oy = pu 3) 
Oy/ox = —(1 + xy)pV (4) 
O¢g/Ox = (1 + xy)u (5) 
Og/Oy = V 6) 


It was shown by Friedrichs® * and by Lighthill® that 
it is advantageous to interchange dependent and inde- 
pendent variables. This transformation from ¢(x, ¥ 
y(x, vy), tox(¢g, v), v(¢, W) is, of course, subject to the 
restriction that the Jacobian of the transformation 1s 
nonvanishing—4.e., 


J = 0¢,~)/O(x,v) = p(l + xy)q? ¥ 0 (7) 


The proposed transformation is therefore valid in the 
entire flow field of the nozzle since only stagnation 
points (¢ = 0) and the centers of the curvature (y = 
— 1/x) are excluded. 

Before proceeding with the transformation it will be 
convenient to introduce independent variable (£, 7), 1 
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TWO-DIMENSIONAL 


lieu of ¢ and y, as follows; 


v = p*q*n (8 


o= fax) dx (9) 


s 


where g)(x) is the known (or assumed) velocity distri- 
bution along the reference streamline yo(x, y = 0) = O 
orn (x, Y = O) = O; and p* and q* are the critical 
density and velocity, respectively, corresponding to 
sonic flow conditions. 

Another useful parameter, /(x, y), is proportional to 


4 one-dimensional area ratio and is defined as 
h= (p*q*) (pg) = A/A* (10 


Substitution of Eqs. (8) through (10) into the trans- 
formed Eqs. (3) through (6) leads to the following four 
differential equations for x, y, g, 8, and h(p, qg) as func- 
tions of € and yn, and the known parameters «(x) and 


q(x : 
ax dE = (Ox /d¢) (dy/dé) = (qv'q) {cos 8/(1 + xy)] 
(11 
Ox/On = —h[sin 6/(1 + xy)] (12 
Oy OF = (qe qg) sin 6 (13) 
Oy/On = h cos 0 (14 


Elimination of x and y (except for the 1 + xy) leads 
to two first-order partial differential equations that are 


nonlinear as follows 


0 On(gn g) + (go'g) (1 + xy) cos? 0(0/0n) [1/1 4 
xy) |] = h\ (1 + xy) sin @ cos @ (0/0) X 
(1 (1 + «v)] + (00/08)f (15 


gv q) ,;(08/On) — (1 + xy) sin @ cos 6(0/ On) 
[1/(1 + «y)]f = (O8/0E) + 
h(1 + xy) sin? (0 0) [1/11 + xy)] (16 
These are, of course, insufficient to determine the 
four unknowns g, y, 8, and #. The additional relations 


needed are Bernoulli's equation 


(q°/2Z T ( 7 l 
Ky + 1)/2(7 — Iie™* (17 
where « V yp p is the local speed of sound, and the 


isentropic equation of state 
Pp = p (1S 


After eliminating the pressure from the speed of sound 
one can write Bernoulli’s equation in the form 
p p™ rly + 1)/2] - 

91 
l(y — 1)/2] (q/q*)?} (19) 
This equation together with Eq. (10) gives the needed 
relation between g and h, 
l 


l(y -— 1 2] (q q*)?t (20) 


Hence, 
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h/hy = (q/qo)' VI(y + 1) — (y — 1) (qo/q*)? X 
(¢ qo)? | l(y + 1) — (y —_ 1) (qo/q*)*]} : (?] 


It is also desirable to determine the Mach Numbers 
M(x) along the reference streamline in terms of ho(x) 


or q(x). From Eqs. (17) and (20) one obtains 
My? = (qo/c)? = [2(qo g*)?|/[(y + 1) - 

(y — | (q q* (22 
hy = (1/My) {2 + (y — Mey + DP OHP/BO- 


For direct comparison with experimental data, the 
equation for static pressure will be useful. The ratio 
of the pressure p(é,n) to the pressure p(é,n Q) is ob 
tained from Eqs. (10), (18), and (19) and is 


* 


P po = [(e'p*) (9/9*) (9* QI" (P*, po) = 
(g q (hih 24 
where py) = P»(x) is the static pressure distribution along 


the reference streamline 


Solution of Differential Equations 


Following the example of Friedrichs® ‘and Nilson,’ 
a solution of Eqs. (12), (14), (15), and (16) is obtained 
by assuming a power series in y for the dependent vari 
ables xv, y, g, and @ and by utilizing the relation between 
g and h given in Eq. (20). 
0 be the wall of the 
nozzle along which the velocity distribution g» and 


Let the reference streamline 7 


curvature «x are known, and let the origin x 0 be 
yy = q* 
0 can be written 


taken at the point where Then the boundary 


conditions along the wall n 


n= 0 
\ 0 a 
23 
f 0 
q q 
Since (kv),.9 = O, the substitution of Eq. (25) into 
Eqs. (11) and (12) gives 
(Ox Of), 0 l | : 
(26 


(Ox On), —o Of 


With the above boundary conditions, series solu 


tions for x, y, g, and @ can be written as follows 


\ E+ x,(E)n + X2(E)n- + x3(E)n°* 4 (2% 
) vi(E)n + WolE)n~ +- y3(E)n” 4 »S 
Gg = ql + qil&)n + g2(&)y- + gsl&)n’ 4 cad (29 
A = O,(E)n + O2(E)n> + O3(E)n? +... (30 


The coefficients of powers of 7 can be found by substi 
tution into Eqs. (12), (14), (15), and (16). The param 
eter / in these equations is written in terms of an 7 
series by means of Eqs. (21) and (29). 

For nozzle applications the curvature is small enough 
so that «x? and higher order terms can be neglected. 
After frequent application of the binomial theorem 
and equating like powers of n, one can solve for x-, y-, 
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g-, and 6-coefficients in terms of the known x(é), /o(&), 
and ./,(&) and their derivatives with respect to & which 
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are denoted by primes. Substitution of these coeff. 


cients into Eqs. (21), (24), and (27) through (30) gives: 


bre! yg {shat Me ape oy g RM? = 1) 
=§- F : Ly“ Sie — 21 4 ; 3] 
Pinot . 7 16(Mu? — 1) °°" ™ 65 fo 
vy = hon — } [who?(My? — 1)]/2}n? + (0/6) [Rohe (Me? 1) hy’? |n 39 
tide ~~ 12 ly(y + DMS — (272 — 1)M 3yMo! — 3M? + 1 
= 1 — kh a vy + 1)Mé — (27? — 1)Mé — 3yMi! — 3M? + 1] + 
qo — 2 " V6(My? — 1)8 °°" : ' 
khy*he’’ y+ 4) Mot V2 ») 4 K'hy*hy’ mie Kh i i)! 
iS eo ae 21 T Ly 7 1) = ’ 3 
6(My?— 1) °° 3(Mie2 — 1) °° " 6 (" 
: Khoho'(y + 1)Mot | w’he?(? — 1) |, fho*ho’’’( Mi? — 1) hho’ hy!’ 
0 = Io'n — 24. 
hon | 2(M.2 — 1) 2 |, v4 6 v6 
(y + 1)1M' kk’ hh / 
| (M2 — 1) T 9 [ l (7 } 
h | F ( Vv és | 4 hohe’ (M0? = 1) ~ j khohy’? 41 | Vl ; Vl 
= = ine _ = a (4 - hd a — (D+y- y . 
ho a mn 2 T YoU? — 12°" ; 
. 2k hy?ho’’ 5 K hy7hy’ KN? j 
37M —_ 3M? + 1] ob [(y + l )Mo4 = M,? + 1] + ‘a (¥ + ) VU, + 6 ( | l n 5) 
” ) 
p Noh" seer ~§ whole’? : 
= 1 + khyMo'n — y Mi? 1? + yMo? - = — [y(y + 1) M8 — (2y? — LANE — By 1h 5 Wy? + 
Po Z (6(M.? = 1)? 
ii ie + ORS a ~ te yates +" Can? — 1)! 
5 bilo KM — 4 ).ALy* FF (by Y 36 
6(My? — 1) ‘" : 3(My? — 1)" * ( j= | 


The above set of equations give, for constant 7, 
the streamlines (x, y), the speed g, the flow direction 
6, the parameter /, and hence the Mach Number J, 
and the static pressure p in regions I and II of Fig. 1. 
All of these quantities are given in terms of the curva- 
ture «x and the /y or J) distribution along the reference 
streamline » = 0. The geometry of the flow field is 
summarized in Fig. 3. 


Discussion of Solution 


With regard to the validity and uniqueness of the 
above solution, one can only conjecture that under cer- 








Variables and parameters of curved nozzle. 


Fic. 3. 


tain conditions a unique solution is likely. Such rea 
soning has been given by Friedrichs’ for the case of an 
axisymmetric nozzle and by Lighthill® for a two-dimen- 


sional nozzle. Their arguments are as follows 


Instead of solving the boundary value problem for 
5) to (6) with prescribed nozzle contours and 


Eqs. 
stagnation conditions, one reverses the problem and 
determines possible nozzle contours from a_ velocity 
distribution go, prescribed on one wall. As long as 
this gq or fo-distribution is analytic, a unique solution 
of the differential equations exists, and the series ex- 
pansion converges in some neighborhood of the wall. 
One, however, does not know the extent of this neigh 
borhood, except that as long as qg or / is smooth and 
analytic along every streamline within the field, the 
chances of having a valid solution are good. Conse- 
quently, one may trust this solution as long as the 
streamlines are smooth analytic curves. If analytic 
streamlines are to be pierced together, one must re 
quire at least continuity in the curvature of these con- 
tours. 

For practical applications one must only watch the 
value of », chosen for the extreme streamline. The 
results of two-dimensional symmetric nozzle calcula- 
tions by this method! indicate that for design Mach 
Numbers less than 3, and with series which are carried 
out to n* terms, 7, should not exceed 0.20 and should 
preferably be around 0.10. The higher the design 
Mach Number, the smaller should be the maximum 
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value of 7, otherwise terms in higher powers of 9 should 
be retained. 

If one calculates a nozzle contour by this method 
ind finds that the 4,-streamline is a smooth and con- 
tinuous curve and that q or / along it are also continu- 
smooth, the validity of this solution to the 


ous and 


po dq hihi ed 
1 + who( Mo? — 1)q — — (Mi? — 2)n? + 
[oo + (472 — By + 3) Mio + Cy Mo! + 4M? — 1] 
14.\/ + ) | + Kho*ho (y + 
3(M,? — 1 


For specific nozzle problems for which the design 
Mach Number Wp > 3 and n, > 0.20, it is suggested 
that Eq. (37) be calculated for n, and plotted vs. & to 
detect any strong trend of the Jacobian toward zero. 
Since the Jacobian applies only to the transformation 
process and not to the series solution of the differential 
equations, this would appear to be a weaker test for 
the validity of the solution than the smoothness test 


of n. and h. 


Special Boundary Conditions 


All of the final equations given in the preceding sec- 
tions contain the parameter /y». It will therefore be 
useful to write the boundary conditions of the nozzle 
problein in terms of fy. From the definition of / it 


follows that 
hy 2 1 


a 


The corresponding values of the Mach Number are 


found from Eq. (22) and physical considerations 


VWn)e=0 °| 
Mo’)g=9 > 0 (39 
Wo) p : u,| 


Since from Eq. (23 


Pio Mo Mo! (My? — 1)]/} Me? [2 + (vy — 1) Me? ]j 
the above conditions imply that 
[ho’ “(M |) |e0 = finite 
hence 
hy’)z.9 = O (40 


Furthermore, all terms in the equation for the speed 
ratio, Eq. (33), must remain finite as —~ 0. The first 
two terms of the coefficient of n* in this equation require 


together with jy)’ = O, that 


specific problem seems to be assured. 


As discussed previously, the nonvanishing of the 
J 


Jacobian of the transformation from (¢, W)- to (x, y 
variables is a measure of the validity of this transfor- 


mation. From Eqs. (7), (10), (32), (33), and (35), the 
Jacobian becomes 
hiyhy’? , 
; = ly(y + 1)M ty 27 1) X 
6(M.2 -— 1 
se [2(2y — 1I).MS — 5(y — 2)Myt — 
6 (M,? — 1) 
; ; ; Kh ai a 
1) (My? — 2)Mo* + (VM 2) (Mo? - 1) en? (37 
6 
[x (VM 1) Jeo finite 
and 
[khy’’ /(.M 1) Jeo finite 
12, 


This last result is interesting in that zero curvature 
at the throat is required without any assumptions hav- 
ing been made regarding the shape of the sonic line; it 
is a consequence of the general throat boundary cond1- 
tion applied .to the n*-terms of the g go equation. It 
should be recalled that Gértler'* has shown that vanish 
ing wall curvature implies a straight sonic line and that 
Sauer!’ has shown that finite curvature at the throat 
results in a curved (parabolic) sonic line. 

An asymmetric curved nozzle with a curved sonic 
line is likely to introduce shock waves right at the 
throat. This condition exists whenever expansion 
waves hit the sonic line where they are reflected as 
shock waves. The curvature of the sonic line is 
rather strong near the expansion wall, and the first 
few expansion waves downstream from the throat are 
likely to hit the sonic line. Thus, it appears that for 
curved nozzles a straight sonic line is to be preferred, 
which is assured if the wall curvature « goes to zero 


at the throat. 


Equation of Patching Mach Line 


The right running Mach line that intersects the start- 
ing contour (reference streamline) at a point where the 
design Mach Number is reached is called the patching 
Mach line. It divides the region of general expansive 
flow from the simple wave flow (region III, Fig. 1 
The latter, in turn, is terminated by the design Mach 
wave where the desired uniform flow is obtained. In 
order to calculate the upper contour, which closes the 
region of simple wave flow, the Mach Number and 
flow direction along the patching Mach line must be 
known. 

The equation of the patching Mach line in (£, 7) co- 
ordinates is determined from the fundamental property 
that the angle between flow direction and the tangent to 


a characteristic (Mach) line at point P is the Mach 
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angle a at that point. From the geometry of Fig. 4 a = sin—! (1/M) 
it is seen that 


B = (4/2) — a = tan! VM? — | 
dy/{(1 + xy)dx] = tan (0 — a) = 


—{(1 + tan @ tan 8)/(tan 6 + tan B)] (41) In terms of the (,) coordinates and the appro. 


priate n-series expansions, Eq. (41) reads 


where 
dn —| ( khy hyho’’ i 
a | (y + 1)Moin — (y + 1)Moin? + (y + 1)M' X 
dé hoV My? — | ¥ 2 (M,? — 1) 1(.My? — 1) : . ‘ 
J hohe * [y(y + IM — (2y? — 1) Mie — 8yMot — 3M? + 1) 4+ hott 
(12(My2 — 1)! 24(My2 — 1)? 
, hy*hy’ ' x" hy?) 
(Sy — 1)My' — 4(3y — 4)M,? — 22 : (+¥M,* + 2M.? — 1) 4 9 
i ' I+ coe — pe Ot T 1 (7) ™ 


The solution of this equation is accomplished by any method of numerical integration. It is stressed, however, 
that this integration should be carried out carefully since the accuracy of locating the reeompression contour will de- 
pend on an accurate knowledge of the patching Mach line. 

The Runge-Kutta method of integration is considered most suitable for this purpose and is briefly outlined below 
Frost!! calculated a patching Mach line for a two-dimensional symmetric nozzle by Runge-Kutta and by character- 
istics and found virtually no error. The series approximation, also given below, however, showed noticeable devi- 
ations from the Runge-Kutta patching Mach line even for a design value of n, = 0.10. 


Runge-Kutta Method 


Equation (42) for the patching Mach line can be written in the form 


dn = [|Ao(é) + AilE)n + Ao(E)y? + A3z(E)n*|déE = FlE,n)dé (43 
where 
Ao = —hp-"(M,? — 1)7" 
A, = —Ut[k(y + 1)Mi4)/ (20M? — 1)77)]} 
Ay = [hy (y + 1) Mi']/ (4M? — 1937] 
(y+ 1)Met (  xho”? , 
, [yy + 1I)MS — 27? — IMS — 8yMot — 382 + 10 + 
’ 6(M,2 — 1)'/? (20.2 — 1)! 
Khyhy [Sy — 1) Me — 4(87 — 4)My? — 22) + — hoho (yMot + 2M? —1)+— hi? \ 
1(My? — 1)? (My? — 1)? ‘4 





The first step in the Runge-Kutta integration starts 
at the known point Po(&,y) = 0) with a suitably chosen 
€-interval Aé and is determined from 


An = (1 6) (A, + 2K» + 2K3 + K;) (44 
where 
ky = F (£0, ) AE 


K» 


II 


Fé + (AE 2), nm + (K, 2) | Ag 
Kz; = Fl(& + (AE/2), m + (Ke/2) AE 
K, = Fé + Aé, No + K3) Aé 


The first point on the patching Mach line, Pi(£1, ), 1 
then known, since 








b= & + Ak 
x m = m0 + An 
oO — 
Ihe same equations apply to successive intervals un 
Fic. 4. Patching Mach line. til » = 7, is reached. 
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Series Approximation 


For small », say for less than 0.10, a fair approxima 
tion for the patching Mach line is obtained by ex 
nanding £ into an y-power series about point P, 


& = do + Qin + Gon” + sn” (45 


The a-coefficients are determined by equating corre- 


sponding -terms of Eq. (42) and the differentiated 
Eq. (49). The results are 
a = (&)p, = §p 
a, = —Inyp! M p? = 1)! . 
h ph p(y . a l ) Mp? ?2( Mp? — | 
dy = af 
: 1( Mp? — 1 
K hy p(y + 1) Mp? 
1(Mp? — 1)'* 
hy plop’ . 
a = => [(y? + y)Mpi — 2y?Mpt — 
] 12(Mp? — 1)? L\'Y Y D Y I 
' ho n*hop”’ 
(7 — 1)Mp* — 6Mp* + 2] - "2". x 
12(.Mp? — 1)°* 
[(y + 1) Mp! — (Mp? — 1)] + 
ho pho p’ ( Mp* 
Khop*hop'(y + 1AM (Sy + 3)Mp* — 
24( Mp? - 1)? 
2(27 4 1) Mp? _ 1] 4 K hop*(y + 1) Mp? 
12(.VUp? — 1) 


Once the coefficients a, a,, a2, and a3; have been deter- 
mined, the parameter & for the patching Mach line can 
be calculated for any value of 7 from Eq. (45). 


Calculation of Simple Wave Region 


Once the patching Mach line my is calculated, the 
flow inclination and Mach Number along it can be 
found. These quantities are sufficient to determine the 
nozzle recompression contour for the simple wave 
flow region, which turns the flow from patching Mach 
line m, into the desired uniform test section flow start 
ing at the design Mach wave m. 

From Fig. 5 and the continuity equation expressed in 


stream function terms one gets 
Yo -— ¥=Vr— v= pel 


where c is the speed of sound at P, in a direction, by 
definition, normal to the characteristic line PR. 
The length of the characteristic line is 


l= [(¥o — ¥)/pc] = [(e*q*)/(pc)] (ng — 4%) = 


(hM)p(ng — 9) (46) 


Since h = A/A* 
tropic table,'* it can be read if either VM, M*, p/H or 
any other flow property of point P is known. 

The length / of the characteristic from P determines 
point R on the nozzle contour 


is tabulated in any supersonic isen- 


i.€., 


R= R(X + /cos(at+60+ wo), V+ 
lsin (a + 6+ w)] 


(47) 


where a is the Mach angle at P. 


ASYMMETRIC NOZZLES 








Recompression region 


Fic. 5 


In this manner as many points as desired may be 
calculated for the recompression contour between Q 
and S. 


Initial Mach Number Distribution and Wall Curvature at 
the Patching Mach Line 

The equations derived so far give an insight into the 
relation between design Mach Number, its gradient, 
and the curvature of the contours at the design Mach 
wave and the patching Mach line. 

At, or near, the point P), where the initial Mach 
Number reaches its design value, the curvilinear co- 
ordinates x and y can be approximated by Cartesian 
coordinates « = & 

The differential expression for the curvature of a 


— and ¥ = n (Fig. 4 


streamline 7 = F(#) near point P, is therefore 
kp, = lim } (d2y (dé) [1 + (dy dé)? ]" “ 
nO 
which can be written from Eq. (45) as 
Kp holo’ ([(y + 1) Mot — 2 (Mi? — 1)) 
P21 + Ay? My? — 1 eb 1S 
The wall curvature is zero downstream of Po and 
smoothly approaches zero from upstream, thus 
» =O0 
hy m 
and 
1» = 0 
Mo, 
all other terms of Eq. (48) being finite at P;. Thus, 


the Mach Number distribution ./), initially assumed 
along the lower contour, must have a continuous gra- 
dient that goes to zero at the patching Mach line—. e., 
(dM /dx)p, or (dM dé)p, = 0. The same statement, 
of course, applies to the distribution of the parameter 
hy and becomes an important boundary condition for 
the fy = ho (£) distribution. 

The fact that the Mach Number gradient at P, 
vanishes assures also continuity of curvature of the 
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Fic. 6. Relation between nozzle translation and design Mach 
Numbers 


opposite contour at the patching Mach line (point Q 
in Fig. 4) and at the test section juncture (point S in 
Fig. 5) (provided of course that dn/dé is continuous). 
This was shown by Frost® and Evvard and Marcus." 


DESIGN MAcH NUMBER CHANGE THROUGH CONTOUR 
TRANSLATION 


It has been shown experimentally that an asymme- 
tric curved nozzle can be used to vary the design Mach 
Number by a simple translation of one contour with 
respect to the other. The design of such nozzles has 
been by the method of characteristics and a rather 
lengthy process of iteration. The question arises: 
Is it possible to derive relations for the analytic nozzle 
design method to include the Mach Number variation 
by means of translation or is an iteration method indi- 
cated ? 

Suppose one has to design a nozzle for a range of de- 
sign Mach Numbers from Wp, to Mp,. The relative 
translation is easily found from Fig. 6 and is 


d = t(cot ag — cot a) = 


t((V Mp? —1— VMp2— 1) (49) 


No difficulty is experienced in determining a nozzle 
for, say, Mp, after choosing a convenient and con- 
tinuous /»-distribution and calculating an upper con- 
tour for a constant nc, Now consider this nc, contour 
translated through a distance d to give, supposedly, 
a design Mach Number MVp,. The 7-value of this 
translated contour is from Eqs. (8) and (10). 


Wc. = ft hoy (50) 


This ne,-streamline contour is supposed to give a 
uniform Mp,-flow. In order to verify this, one would 
have to solve the boundary value problem for Eqs. 
(3) to (6) with prescribed streamline nce, and given 
stagnation conditions. This is exactly the problem 
that has not been solved yet and which is avoided by 
the preceding method. Thus one is unable to deter- 
mine the flow field under these conditions directly, and 
a process of iteration is required. 

Suppose one starts with an assumed /p,-distribution 
of the same form as /y, but with a design value of Roy 


corresponding to V/p,. After calculating an upper con- 
tour with this /4o, function for an nc, value, one finds 
that, in general, this contour will not exactly coincide 
with the nc, contour translated through a distance ¢. 
The test section height is of course correct and gp 
is the throat; then, if the two /»-distributions are cop. 
tinuous and have zero slope at the throat and at the 
start of the test section, continuity of curvature is also 
assured at the patching Mach lines and the test section 
junctures of these contours which are smooth and con- 
tinuous throughout. The discrepancies between these 
contours cannot therefore be very large and can be 
approximated by a perturbation term to the /y,-distri- 
bution. 

Assume the /y,-distribution perturbed by a small term 
6(£) which is continuous and has small, continuous 
derivatives, 6’ and 6’. Substitute 


hy = ho aa 6 
hy’ = In,’ + 6’ 
hy’ = hy!" 1 5!’ 


into Eqs. (30) through (36) and neglect products of 6 
and its derivatives with 7”, n°, and «x. It is found that 
the perturbation terms remain only in the y and 4 


equations as follows: 


v = v(ho,) + 6n 

6 = A hi,) + 6’ (51 
where (/y,) denotes the value due to the unperturbed 
hy,-distribution. 

It is therefore possible to determine a perturbation 
function 6(£) from the geometric discrepancies Ay and 
A@ between the n-.-contour calculated from an /y,- 
function and that of the translated n¢,-contour. A 
second yc,-contour can now be calculated from the 
(fo, + 6)-distribution and compared to the trans- 
lated nc,-contour, Fig. 7. This process may be re- 
peated until agreement is satisfactory. 

Since Ay and A@, and therefore 6 and 6’, will be small, 
this iteration will converge in a few cycles. If the 
method of characteristics is used in the design of such 
a variable Mach Number nozzle, about 10 or more 
iterations are required. * 

It should be noted that Eq. (51) is valid only in the 
nozzle region II (Fig. 1)—i.e., between the throat and 
the patching Mach line. The recompression contour 
(region III) is uniquely determined only by the geom- 
etry of the patching Mach line and the flow properties 
along it and cannot be related explicitly to the contour 
of region II. This dilemma is greatly alleviated if the 
nozzle incorporates continuous wall curvature at points 
QO and S (Fig. 6)—i.e., if top’ and hep’ + 6’ vanish at 
point Py. One is then assured of having continuous 
curvature and slopes at the end points Q and S of the 
recompression region. Furthermore, point S is geo- 
metrically fixed by the design Mach wave and the test 
section height, point Q is approximately located 
through Eq. (45) and the recompression contour is a 
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Although theo- 


an infinite number of recompression contours 


snooth curve through these points. 
retically 
could exist, the above conditions considerably limit the 
number of curves suitable for a variable Mach Number 
nozzle, and only a few iterations of the type described 
above will be needed to arrive at the proper contour. 
5 
BOUNDARY-LAYER CORRECTIONS 


Boundary-layer corrections to the perfect fluid con- 
tours may be applied to this type of nozzle using the 
conventional techniques.” The nonviscous nozzle 
contours are generally ‘‘opened up” by an amount equal 
to the displacement thickness so that the resulting test 
section flow has the exact design Mach Number and 
is of uniform quality. In the case of an asymmetric 
variable Mach Number nozzle, however, the need for 
achieving the exact design Mach Number no longer 
exists since an additional translation will increase the 
actual Mach Number to the desired value. Poundary- 
layer corrections may still be required to avoid non- 
uniformities in the flow due to viscous effect. Until 
the magnitude of the viscous effects has been determined 
experimentally, it is suggested that asymmetric Mach 
Number nozzles with a wide range of Mach Numbers 
be corrected for boundary-layer growth. 
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Iteration in Semidefinite Eigenvalue Problems 


B. M. FRAEYS pe VEUBEKE* 


Unwersiy of Liege, Belgium 


SUMMARY 


General rules are derived for setting up iteration matrices yield 
ing convergence towards the elastic modes of semidefinite self 
adjoint vibration problems. They are classified according to the 
methods used for a preliminary reduction of the problem to a 
finite number of degrees of freedom. These include the Rayleigh 
Ritz, the complementary energy, and the integral equation ap- 
In the last case it is shown how extended influence co- 
A numerical 


proach 
efficients are related to isostatic reference frames 
example is treated and a procedure suggested for improving the 
higher frequency modes obtained from the lumped mass models 
which are provided by the application of numerical integration 


and collocation to the integral equation 


INTRODUCTION 


— VIBRATION modes of continuous elastic 
structures form an infinite set with a discrete fre- 
quency spectrum. Restraining the freedom of the sys- 
tem to a few of the normal coordinates offers an ade- 
quate approximation to its elastic and inertial behavior 
in some frequency range. Most of the requirements 
of aeroelasticity are for the low range from zero to some 
limit above which the higher frequency modes may be 
assumed not to be appreciably excited. 

It is always possible in theory to obtain the modes 
and their eigenvalues from the differential equations and 
boundary conditions of the problem. This approach, 
however, is impractical except in the simplest situations 
where the general solution of the differential equations 
In this connection a new method due to 
As an alternative, 


are known. 
W. E. Milne! deserves attention. 
iteration methods are useful but present difficulties in 
the case of unrestrained systems, like structures of air- 
craft in free flight. 

Such structures have no Green's function or static 
influence function in the accepted meaning of displace- 
ment at one point due to unit load applied at another, 
and their vibration problem lacks the usual integral 
equation formulation. 

To each degree of freedom remaining in the structure 
after complete rigidification, there corresponds a vibra- 
tion mode of zero frequency. More generally, the pres- 
ence of eigensolutions with zero eigenvalues, which 
characterizes semidefinite eigenvalue problems, _re- 
quires an extension in the concept of the Green’s func- 
tion. 

From this standpoint, integral equations involving 
extended Green’s functions are briefly treated by 
Courant and Hilbert; an example in the aeronautical 
field is given by R. L. Bisplinghoff, G. Isakson and T. H. 
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H. Pian.* H. A. Fettis' applies iteration directly to the 
differential system of those problems. 

Setting aside the very exceptional cases where the 
integrations may be carried out in closed analytical 
form, both of these approaches imply some type of 
numerical integration, whereby a system of linear alge- 
braic equations approximates either the differential 
system or the integral equation. In the present paper, 
reduction to a vectorial space of finite dimensions jis 
accepted as the initial step from which developments 
necessary to provide convergent iterations are left to the 
care of matrix algebra. 

This approach allows in particular a general formula 
tion of the ideas on which extended Green's functions 
are based. The various rules obtained remain of stand- 
ard application in very general situations; all have in 
common the central role played by a matrix operator 
for orthogonal projection in the vectorial subspace of 
the elastic modes. 

MODES 


(1) EQuATIONS OF MOTION AND NORMAL 


Consider the eigenvalue problem exemplified by the 
natural modes of vibration of a system having a finite 


number of degrees of freedom 
(C — w#*M)x = 0 


where C stands for a symmetrical semidefinite positive 
matrix of stiffness coefficients, .1/ for a symmetrical 
positive definite matrix of inertia coefficients, w for the 
circular frequency and x for the one column matrix or 
vector of the vibration amplitudes. 

Eq. (1) is obtained from the general equation of mo- 


tion 
Mq +Cq=p (2 


by letiing the vector / of generalized forces, correspond 
ing to the generalized coordinates vector q vanish and 
by separation of the time variable g = ey. 

To each eigenvalue w,”, root of the frequency equation 


det(C — w*M) = 0 (0 


there corresponds a normal mode .x;,) satisfying the rela- 


tion 
Cx = o,?Mxy } 


The assumption that C is semidefinite is equivalent to 
stating that Cg = 0 admits of non-zero solutions. Let 
ui) (2 = 1, 2 .m) denote a complete linearly inde- 


pendent set of these free modes so that 


Cui, = 0 (@=1,2...m) (0 
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From Eq. (4) it appears that the free modes are nor- 
mal modes of zero eigenvalue. They involve no elastic 
deformation energy and are consequently easy to deter- 
mine by simple inspection of the degrees of freedom re- 
maining in the structure after complete rigidification. 
Elementary examples of free modes are the vertical, 
longitudinal, and lateral translations, the pitching, 
rolling, and yawing motions of a rigid airplane in free 
To those might be added the relative control 


flight. 
surface motions provided they are not effectively 
spring-restrained. 

The notation x,,, (r = 1, 2...) will qualify elastic 


modes only; those are modes calling for elastic deforma- 
tion energy and consequently associated with non-zero 
eigenvalues. The following orthogonality relations are 


essential :* 


ui)’ Muy = 0 1 Fj (6a) 
ui’ Mxy Xi)’ Muy = 0 (Gb) 
Xn’ Mx is) = 0 r#s (6c 
Kiss tx =f ¥#S (6d 


They hold as a consequence of Eqs. (4) and (5) between 
any two modes of different eigenvalues and may be as- 
sumed to hold between two distinct modes of the same 
eigenvalue as the result of an orthogonalization process. 


The modes 
u ee]: ?...... 8) Xi 


then form a complete orthogonal set in terms of which 
any vector admits a unique expansion. For example, 
the expansion of the displacement vector 


n 
- > EX (7 
I 


m 


q = YS nite 
1 


introduces the normal coordinates 7;(¢) and &é,(¢) and 
may be substituted in the equations of motion Eq. (2). 
Multiplying then to the left by uw)’ and x;,)’ in turn, 
using Eqs. (4), (5), and (6), the equations of motion are 


obtained in normal form 


Ui) P Ss 
"= : (4a) 
ui)’ Muy 
of . X(r)'P ‘a 
w,-& + & = (7b) 
' i Xr) ' Mx 


The simplification resulting from the absence of inertial 
and elastic coupling between normal coordinates justi- 
fies the interest attached to the computation of the 
normal modes. 


(2) PROJECTION OPERATORS IN RESTRAINED SYSTEMS 


Let the normal mode problem of a restrained system 
be cast in the form 


(C-'M — rXE)x = 0 (8) 


“The exchange between rows and columns or transposition is 
denoted by a prime. wu’ and x’ are therefore one-rowed matrices 
T . 
in transposing a product each factor is transposed and the order 


reversed 
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where C~' is the reciprocal of the nonsingular matrix C, 


/ the unit matrix, and 


d= 1/e? (9 
the characteristic value parameter. With the trans 
formations 

M = LL’ L’x =y 10) 


It reduces to the characteristic value problem 
(F — AE)yvy = 0 


of a symmetrical matrix 


The computation of Z is always possible and elemen- 
tary when it is assumed to be lower (or upper) triangu- 
lar, L’ then being upper (or lower) triangular. This tri 
angular factorization has been proposed by Cholesky in 
connection with the solution of symmetrical linear equa 
tions systems.° 
be the modal columns and A; > 2. 
To eliminate the mode ya 


Let > Ay 
the characteristic values. 
yielded by the iteration with F we form the scalar prod- 
uct Vau)’Vva), the degenerate square matrix ya)Va)’, and 


the resulting projection matrix 
A, = E — [ywya’/yay’ 
In view of the orthogonality relations between modes 
ya Vn = 0 yr #1 
this operator has the obvious selective properties 


JAiya) = 0 


(Ain = Yo rl 


We have accordingly for the deflated matrix® 


F, = AiF = FA, = F ee Aly pa , Ya 'y 1 
the relations 
J Fy 1 = O 
l Fy = Avo ry #1 


proving that F, retains the modes of the original / 
matrix. (r * 1) retain their 
original characteristic values \,, the characteristic value 
The very first iteration 


But, whereas the \, 


of ya) has been reduced to zero. 
with the deflated matrix /; will then remove any first 
mode component and will converge toward the mode of 
next highest characteristic value. Pre- or postmultipli- 
cation by the successive projection matrices allows re- 
petition of the iteration process without the growing 
complications sometimes encountered with other 
methods like ‘“‘sweeping.’” 

The transformation (10) is not an essential step and 
may be avoided altogether. The projection matrix to 


be used is then 


with the selective properties 
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Axa = 0 Atay = x r l m Uni M 
sia A=E-) (i 


resulting from the orthogonality relations Eq. (6). 
The original problem Eq. (S) is replaced by the prob- 


lem 
(C—'MA; — XE)x = 0 


in which \ = 0 is now the characteristic value pertaining 
tox 1)° 
When interpreted as the substitution of a deflated 


inertia matrix 


Mx) [Mx ]’ 


M, = MA, = Ai'M = M 
Xa)'Mx l 


in place of the original 7, the method reduces to zero 
the effective inertia of the mode to be eliminated. 


Va ‘Mx = Xa)’ MA xq) = 0 


thereby raising its frequency to infinity, while the other 
effective inertias are left unaltered 


tm’ Minn = Xq)'MA Ky = X~y!/Mxq) y sey 
(3) INERTIAL DEFLATION OF FREE MODES AND THE 


RAYLEIGH-Ritz METHOD 


The Rayleigh-Ritz method uses a finite expansion in 
assumed displacement functions for the continuous 
structure, each multiplier being conceived as a general- 
ized coordinate. The strain energy is computed from 
the strains obtained by differentiation of the displace- 
ments and is ultimately a quadratic homogeneous form 
in the generalized coordinates with matrix C. The 
kinetic energy is a quadratic homogeneous form of 
matrix .\/ in the time derivatives of the generalized co- 
ordinates. As a consequence the vibration problem is 
obtained in the form presented by Eq. (1). From the 
positive definite nature of the kinetic energy the matrix 
J\/ is never singular and suggests the use of an iteration 


process based on the modified equation 
(M-'C — wk)x = 0 


with inherent convergence toward the mode of highest 
frequency. Several procedures are available® for shift- 
ing convergence to the lower modes. Wielandt’s 
method for instance® is applicable and enjoys the ad- 
vantage of determining the modes independently and 
unaffected by round-off errors due to the elimination of 
other modes. It requires, however, some a_ priori 
knowledge of the eigenvalues and costly matrix inver- 
sions. In the case of unrestrained systems the more 
promising formulation Eq. (8) is unavailable, since, in 
view of Eq. (6), the C matrix is singular and its recipro- 
cal, the matrix of static influence coefficients, does not 
exist. 

However, such a formulation exists in the vectorial 
subspace of the elastic modes and will be shown to follow 
from inertial deflation of the free modes. 


The projection matrix used here is 


T my’ Mu 
with the selective properties resulting from Eqs. (6 
Au, = 0 (a@=1,2...m) (12 
AX) = XG (ee 1.2)... 8) (13 


In view of Eqs. (13) the elastic modes of the modified 


eigenvalue problem 
(C — w2°MA)x = 0 (14 


are identical in shape and eigenvalue to those of the 
original Eq. (1). On the other hand, the free modes are 
now modal vectors of Eq. (14) with arbitrary eigen- 
values, since they satisfy both Eqs. (5) and (12 

There follows that .1/A presents the same type of de- 
generacy as C, which indicates the possibility of a simul- 
taneous reduction to the vectorial subspace of the elastic 
modes. For definiteness we assume that in the applica- 
tion of the Rayleigh-Ritz expansion to the unrestrained 
structure a complete set of m independent rigid body 
displacements is used as representing the m first as- 
sumed functions, the 7 remaining ones involving com- 
patible elastic deformations. 

The free modes are then represented in the vectorial 
space by the unit vectors e;;, (7 = 1,2...) defined 
as having all their elements equal to zero except the ith 
one which is equal to unity. Any orthogonalized set of 
free modes, from which the projection operator A 
should happen to be constructed, is made of linear 
homogeneous combinations of these unit vectors with a 
non-zero determinant. Hence from Eqs. (5) and (12 
follow the equations 
Ce, = 0 Aen = 0 MAe.; =0 t= 29. 
indicating that the m first columns of the matrices C, A 
and ./A are composed of zeros. Since further C and 
MA are symmetrical the matrices are partitioned as 
follows 


‘ 00 - OR on 00 
OK “QB ~ — lgln 


K and N are symmetrical nonsingular (”, m) matrices, 
R and B respectively (m, ) and (7, m) matrices. 
It is then natural to consider a similar partitioning 0! 


the modal vector 


where v contains the m first elements, z the m remaining 
ones. Eq. (14) thereby reduces to the elastic subspace 
formulation 

(K — w*N)z = 0 
identical with that of a restrained system. Iteration 
may now be performed with A~“N to obtain the lowest 
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frequency elastic mode and successive deflations applied 


on the basis of the orthogonality relations 
‘N23, = O Zr) Ke 0 rts 


holding in the subspace. The m first components of a 
modal vector may be obtained from the orthogonality 


onditions 


é Mx Q ] 


SS eee 


implied by Eqs. (6b), and also more directly from the 
selective property expressed by Eq. (15) in the form 


Rz ls) 


Ricip Bopy EQUILIBRIUM CONDITIONS AND THI 
COMPLEMENTARY ENERGY AIETHOD 


It is a feature of the complementary energy method 
that internal and external equilibrium should be verified 
. priori. No difficulties are to be expected from unre 
strained systems since the rigid body equilibrium con- 
ditions, implied by each free mode, are used to eliminate 
the free mode coordinates. It is however instructive to 
show how the use of the projection matrix A is effective 
in carrying out the elimination. 


The static problem 
Cq = p (16) 


to which Eq. (2) reduces on the assumption of zero 
accelerations admits of a finite solution for g, but under 
restrictive conditions. From a classical theorem of 
algebra" it is necessary and sufficient that p should be 
rthogonal to all the solutions of the homogeneous ad- 


joint equation 


Since in this case C = C’ the conditions are explicitly 
tan Pp = p'uy = 0 ce L2....m) GHA 


hey state that the virtual work of applied forces must 
vanish when the displacements are in the ratios given 
by a free mode and they are consequently the expres- 
sions of all the rigid body equilibrium conditions of the 
problem. Now, since .\/ is a nonsingular matrix, the 
vectors Ju; and .Wx,,, each one proportional to the 
inertia force distribution associated with the mode, are 
linearly independent and form a complete set in terms 
of which any arbitrary vector p admits of a unique ex- 


pansion 


b= YD aMuy + dX B,Mxy (18) 


Premultiplying by u;,)’ or x)’, and using Eqs. (6) 


again, the coefficients of the expansion are found to be 


Uy) 'p Xin’ P 
a;= ‘ Bb, = : (19) 
ui’ Muy Xin’ Mxy, 


When these values are substituted in Eq. (18), and since 
P 1s arbitrary, the following useful expansion is ob- 
tained: 


\/ , , 

5 Mui yu ~ Mxin3 

I SN : i > 0) 
") fray Ve — Wa 

Eqs. (19) show the conditions (17) for static equilibrium 

are equivalent to 


a 8) ee ml 


They require that p should contain no inertial compo 
nent from the free modes. Hence, starting from an arbi 


trary p, the vector 


p> —dYiaMu (7: z a +) 


u Wu, / 


is the most general one fulfilling the requirements 
This result is susceptible for a well-known dynamic in 
terpretation.” On the assumption that the system 1s 
completely rigidified, the displacements are restricted to 
free modes and the equations of motion to Eqs. (7a 
Now, from Eq. (21 


A'p=p-) ——— Mu 


and in view of Eqs. (7a 


A’p = p a Vu 
The fact that A’p satisfies the rigid body equilibrium 
conditions then clearly results from an application of 
d’Alembert’s principle to the rigidified structure; or, to 
borrow a statement from reference 3, the application of 
A’p is equivalent to ‘‘consider the unrestrained system 
to be in equilibrium such that the external load is 
balanced by the inertia forces corresponding to the 
acceleration of the system as a rigid bodv due to the 
applied load.”’ This new property of the matrix A may 
be used to advantage when setting up the vibration 
problem with the complementary energy method."! 
The same displacement functions are assumed as for the 
Rayleigh-Ritz method, and the structure is loaded by 
the resulting inertia forces, whose amplitudes for har 
monic motion have the vectorial representation 


bp = w*Mx 
We multiply to the left by A’ 
A'p = wA'Mx = w*M(Ax 


and conclude that, whatever be x, the amplitude vector 
Ax is such that the resulting inertia force distribution is 
in equilibrium. Hence taking for x in succession all the 
unit vectors, what amounts to take for Av the successive 
columns of A, we explore all the possibilities for equi- 
librium. Now, from its construction, A has exactly the 
degeneracy m and its m non-zero columns yield exactly 
n independent inertia force distributions satisfying rigid 
body equilibrium. For each of those the corresponding 
stress distribution should be computed by integration of 
the internal equilibrium and compatibility equations. 
From the stresses the maximum strain energy follows, 
whose expression for an arbitrary initial amplitude 
vector x will be of the form 
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(1/2)w*x’C*x 


The element common to the 7th row and jth column of 
the elastance matrix C* is e;;,'C*e,;,, and belongs by 
definition to the interference energy (or energy itself if 
7 = j) between the stress distributions due to the dis- 
placement amplitudes Ae,;, and Ae. Observing that 
the first m stress distributions vanish identically to- 


gether with the m first columns of A 
Ae, = 0 pe 1.9. OO 
the first m rows and columns of C* must also vanish and 
0 0 


C* = a 
OA* 


where A* is a symmetrical, nonsingular (7, 2) matrix. 
For the maximum kinetic energy we have 


2)w*x’A'MAx = 
(1/2)w2x’ MA2x 


(1/2)w?(Ax)'’M(Ax) = (1 


However, as a result of its definition and selective 
properties, a projection matrix enjoys the property of 
idempotence A? = A, whence the maximum kinetic 


energy is also (1 /2)w*v’ MAx. 


We apply now that form of Hamilton's principle 
which asserts that the difference between maximum 
strain energy and maximum kinetic energy is stationary 
with respect to arbitrary variations in the displacement 
amplitudes and obtain the vibration problem in the 


form 
(w*C* — MA)x = 0 


In view of the structures of the matrices involved the 
problem is again carried over in the subspace of the elas- 
tic modes: 


(wA* — N)z = 0 


In order to obtain the low frequency modes, iteration 
should be conducted with the matrix V~'A* and de- 
flation set up with the aid of the elastance matrix A™*. 
It is known from general theory that the complemen- 
tary energy method yields better eigenvalues than the 
corresponding Rayleigh-Ritz procedure.” 


(5) EXTENDED INFLUENCE COEFFICIENTS 


From the considerations developed in section (4) it 
appears that the static problem 


Cq = A'p (22) 
has always a solution for arbitrary p. The general solu- 
tion is 


m 


q = Pop + Sym (23) 
1 


where the first term represents any particular solution. 
The arbitrary constants y; may be represented by the 
scalar product of » with arbitrary vectors g; 


Yi = g;'p 
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so that the general solution takes the form 
q = (P, + Vug ‘\p = Pp (24 
l 


P is the general matrix of extended influence coefficients 

The substitution of Eq. (24) in Eq. (22) leads, since » 
is arbitrary, to the following fundamental equation, 
satisfied by extended influence coefficients 


CP = A’ or P'C=A (25 


It should be compared with the equations 
CC'=C'C=E 
holding for restrained systems. 

A new formulation of the vibration problem, extend- 
ing to unrestrained systems the properties of Eq. (8), is 
obtained from Eq. (14) by left-multiplying with P’ and 
using Eq. (25) 

(E — w*P’M)Ax = 0 
or, in a form adapted to iteration 


where y= Ax (26 


(P’M — r\E)y = 0 


In view of Eqs. (13), the elastic modes are retained with 
their original characteristic values. The original free 
modes however degenerate because of Eq. (12) in the 
trivial solution y = 0. It may be suspected, and ex 
amples do confirm, that the place of free modes is 
generally taken by parasitic modes. It is most im- 
portant that the parasitic modes be eliminated in the 
iteration process and this may be achieved by a careful 
choice of the g; vectors in relation to the particular 
matrix Py adopted. For example, we may require that 
a set of m parasitic modes w,, 
(26) with zero characteristic value, and, therefore, ob 


be eigensolutions of Eq. 


tain as conditions 


P’ Mw,j;) = Po’ Mwy + yb [a ;)’. Mwy |g = O 
1 l 
(g=1,2...m™ 


So long as the w,,, and x;,, are linearly independent this 
yields nonsingular linear systems of the mth order for 


the m X n components of the g; vectors. 


(6) THE SYMMETRICAL EXTENDED INFLUENCE 
COEFFICIENTS 
As a particular choice, the parasitic modes of zero 
characteristic value may be the free modes themselves, 
in which case we obtain, on account of Eqs. (6a). 
P,y'’Mu 


ui)’Mu, 


i= 


The resulting matrix of extended influence coefficients is 


” Ucyt’M 
G« (z ne oe ) P, = AP, 
1 


ui)’Mui; 


It must be independent of the particular choice P». 
Indeed, we have the more general result 


AP =AP,=G (27) 
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which follows immediately from the definition of P in- 
cluded in Eq. (24) and from Eqs. (12). 
It will now be proven that G is symmetrical. As a 


particular P matrix it satisfies both Eqs. (25) and (27) 


G'C . A AG G 


ind the refore 
G’CG _ G 
Now the left side is obviously symmetrical, which 
proves the proposition. 
When symmetrical extended influence coefficients are 
employed, Eq. (26) specializes in 
(GM — rXE)y = 0 (28) 


[his is the most natural extension to unrestrained sys- 
tems of the formulation Eq. (S) an the associated itera- 
tion procedures successfully applied to restrained sys- 
tems. Compared with the original Eq. (1), its eigenso- 
lutions and eigenvalues are unchanged except for the 
free modes whose frequency has been shifted from zero 
infinite characteristic value) to infinity (zero character- 
istic value). 

An explicit form of the extended influence coefficients 
may be found in the case of the Rayleigh-Ritz procedure 


described in section (3). Partitioning the vectors 
» 


to match the partitioned forms of the matrices C and A, 


Eq. (22) reduces to 
Kb = R’r + B's 
lhe general solution for g is consequently 
ja arbitrary 
(b = K-'R’r + K-'B’s 
ind a particular extended influence coefficients matrix 
18 
7] 0 
P. 
te) ee 


From Eq. (27) follows then 


RA~'R’ RK~'B’ 
BK R’ 5K 3p" 


(29 


Nothing much is gained in this case. Actually the itera- 
tion with A~'NV is more economical, both as an initial in- 
vestment and because it proceeds in the subspace of 
elastic modes only. The use of extended influence co- 
efficients becomes a practical proposition when it is com- 
bined in the integral equation approach with colloca- 
tion and the use of an isostatic reference frame. 


‘) COLLOCATION AND IsoSTATIC REFERENCE FRAMES 


Imagine any group of constraints applied to structure 
‘0 as to make it isostatically restrained with respect to 
anigid reference frame. The relative displacements of 
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the structure are then given by 


qd G oP 


1 


where Gj,,. 1s a symmetrical matrix of generalized in 
fluence coefficients in the usual sense of relative elonga 
tions in the generalized coordinates due to unit general 
ized loads. We know that if p is replaced by A’p the 
applied load fulfills the requirements of rigid body equi- 
librium. Consequently, the constraining loads with re 
spect to the isostatic reference frame must vanish. 


Hence 


q Giso,A ‘p 


is a possible displacement vector of the unrestrained 
structure. The role played by the reference frame is 
thereby reduced to determine unique values y,, depend- 
ing on the nature of the constraints, for the otherwise 
arbitrary free-mode components of the displacement [cf. 


») 
( 


the general solution Eq. (23)]. It may be concluded 


therefore that 
P Gig A’ (30) 


is a particular matrix of extended influence coefficients 


and 
G AG;,,, A’ (31) 


the obviously symmetrical one. These expressions fur 
nish us with practical means for computing the ex 
tended coefficients. From them we obtain the formula- 


tion 

(AGi,,A'M — XE)y = 0 (32) 
derived from Eq. (2S), and the simpler one 

(AGiso.M — AE)y = 0 (33) 


deriving from Eqs. (26) and (30). 

In the last case the parasitic modes must be shown to 
be harmless. They are, because it happens that their 
characteristic value is zero. To prove this it suffices to 


show that the equation 
AG. Vy = 


possesses m linearly independent solutions for y; be- 
cause of the existence of 7 elastic modes in Eq. (33) it 
cannot have more than m. Now, from a classical 
theorem of algebra,’ the adjoint equation 


MGiso. A’y = 0 

should then also have m independent solutions; but 
this is obviously true since 

A’'y =0 
has m solutions corresponding to the m solutions 1;;) of 
the adjoint equation 

Ay =0 

Examples of lumped mass systems treated in a man- 


ner equivalent to Eq. (33) are given by S. Levy’ and by 


Scanlan and Rosenbaum.!* In the notations of this 


paper their mathematical treatment is along the follow- 
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ing lines. The equation 
m 


x = w*G,,, Mx + z. a jU¢;) (34) 
I 


expresses that the absolute displacement amplitude due 
to the application of the vibratory inertia loads differs 
from the amplitude in the reference frame by a set of 
free-mode components. The intensities a; of these 
components are determined from the conditions 
wu, Mx = 0 (f= 1.2...) 

of rigid body equilibrium for the applied loads. They 
yield 

, UG)’ MGiso,. Muy 


Uy’ Mu) 


== @& 


and, when substituted in Eq. (34) 


m 1 AY, «)'M 
x= w* E = ps | Gu Mx 


T uj) Mu; 


which brings us back to Eq. (33), 

The formulations Eqs. (32) and (33) are specially in- 
teresting when used in conjunction with collocation. 
Whatever be the kernel considered for the integral 
equation of the problem (usual Green’s function in the 
isostatic reference frame or extended Green's function 
for the unrestrained system), any type of numerical 
integration, using the values of the integrand at several 
division points, and collocation at the same points, 
yields an equivalent lumped mass model with a diagonal 
mass matrix. The elements of G;,,. are then local 
values of the usual Green's function in the reference 
frame, whose definition may be chosen in order to ob- 
tain simple checks against experimental data. The 
elements of G;,,.A’ or G are local values of the corre- 
sponding extended Green’s functions. As a final prac- 
tical observation, the invaluable Rayleigh quotient 


_ x Mx 
x’ MGMx 


is seen to be identical in view of Eqs. (31) and (13) and 
of the symmetry of 7A to the simpler one 


x’ Mx x’ Mx ba 
= = (39) 
(Ax)'’MGig, M(Ax)  x'MGig, Mx 
This makes the computation of G unnecessary when Eq. 


(33) is employed for iteration. 


(S) NUMERICAL EXAMPLE 


The example here concerns the normal vibration 
modes of a free elastic beam of constant cross section 
carrying a concentrated mass at the left end, equal to 
the total mass of the beam [see Fig. 1 (a) ]. 

From the differential equation and boundary condi- 
tions the exact frequency equation is found to be 


u(sinh w cos wh — cosh w sin pw) = 1 — cosh uw cos pu 


The mode shapes are given 


in which pt = w*ma*/EI. 
by 
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Fic. 2. Comparison of second elastic mode shapes on the basis 


of same effective mass: 


= exact mode 
@ = lumped mass model 
© = lumped mass model after recomputation of the effective 
mass 


y(¢) = B(cosh pe + cost) — A (sinh we + sin uf 


where ¢ is a dimensionless coordinate running from —] 
at the left to zero at the right, and 


A(sinh » — sin n) = B (cosh wp — cos u 


The effective mass of a mode may be computed from 
70 
[y2(—1) + / y*(o)de|m = nm 
J -1 


with n = B? + (1/4)y?(—1) 


The two first elastic modes have the following charac 
teristics from which the approximate results found be 
low may be evaluated: 


First mode 





B = 4.041835 ut = 266.878 n = 40.96494 

A = 6.290534 B = 6.321783 

normalized to have y(—1) = 2 
Second mode 

nu =7.1340 pt = 2590 n = 29.43225 

—A = 5.402871 —B = 5.402060 

normalized to have y(—1) = 1 


This mode is represented by the continuous line of Fig. 
9) 


(a) Rayleigh-Ritz 
All boundary conditions are for stresses, none !0r 


displacements. The assumed functions expansion 
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y=qn + af + q3s(2 — 3f + £°*) 


features the translation free mode, the pitching free 
mode about the right end, and the elastic displacement 
due to a load applied at the concentrated mass, right 
endclamped. It satisfies but one of the stress boundary 
conditions: y”"(—1) = 0. 

The result is bad both for mode shape and frequency 


parameter 
wi = 575.34 


This serves to indicate that, although not required by 
the variational principle, good results with few assumed 
functions can only be obtained by satisfying also most 
of the stress boundary conditions. From other pub- 
lished results” this is a critical point to watch in all 
cases of unrestrained systems. 

The next expansion tried, 


y= + at + as (F — 2F* +) 


satisfies both boundary conditions for flexural moments 
but not for shear loads. The elastic part is the deflec- 
tion due to uniform loading of the beam on two end sup- 
ports. 

The main results are contained in Table 1. The error 
on the frequency parameter is about 5 per cent. The 
mode shape is here proportional to the last column of A. 
It was first chosen to be normalized like the first exact 


mode 
y = 2 — 38¢ + 100¢% — 50¢4 


It compares well and even slightly better if the com- 
parison is made, as it should be, on the basis of the same 
effective mass. For this purpose it should be altered in 


Fee sat 
= 0.9759 
13.01587 


the ratio 


TABLE 


(b) Complementary Energy 


The inertia loads due to the same expansions as in 
the Rayleigh-Ritz case were also investigated. Despite 
the bad mode shape the first expansion yields a fair 
approximation to the first frequency parameter ut = 
277.9. Here are some more details for the second case. 

The inertia loading, satisfying rigid body equilib- 
rium, is found according to theory by using the last 
column of A in Table | as vector of displacement ampli 
tude. The load consists of a distributed part 


w*(m/a)(—0.04 + 0.76 — 2¢% + ¢4 
and a concentrated part —w*m0.04 at the left end. 
Integration is carried out twice to find the bending 


moments 


38 I 
ma | —0.04¢ — 0.02¢? 4 fm — Ol 4 ;* 
300 50 


and as a check it is verified that the stress boundary 
conditions may be satisfied at both ends. From the 
strain energy 

_,  wim'a® | 86 SST 
K* = 10 


7 - and yu‘ 
| 135 135 


10.84 135 135 Da 
a 207 61 
630 S66 SS7 
which is correct to about 0.27 per cent. Since only one 
elastic degree of freedom was retained, the mode shape 
and its effective mass are here the same as in the Ray 
leigh-Ritz case. 
(c) Collocation and Isostatic Reference Frame 
If one aims at several elastic modes this method is 
attractive on account of its low initial investment in 
numerical labor. Two elastic modes, besides the two 
free ones, were aimed at out of five degrees of freedom. 
The integral in the integral equation was approximated 


l 


Numerical Results for First Elastic Mode from Rayleigh-Ritz Method 
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by Simpson’s rule with four intervals and the subdivi- 
sion points used for collocation. This is equivalent to 
the lumped mass model of Fig. 1 (b); the coordinates 
are the displacements of the concentrated masses. 

Iteration was performed according to Eq. (33) and 
Table 2 gives the main results. The G;,,. matrix was 
computed from the deflections of the beam on two end 
supports. This causes the first and last columns of the 
iteration matrix to vanish and allows the temporary 
discarding of the first and last coordinates during the 
iteration process. Besides initial investment this is 
another advantage of Eq. (33) over Eq. (32). In that 
respect, clamping the beam at one end would have 
proven slightly less advantageous. 

The frequency parameter is correct to about 0.6 per 
cent; the mode shape (of same effective mass) within 
plotting accuracy of the exact first mode. 

The results of a subsequent deflation and iteration 


, 9 O 12 IS 
yh 0 @& 

Fi eT 0 

deen 0 -12 —18 

O -l5 —24 


The application of the formula gives the mode values at 
mid-stations: 
13.013875 — 3.514210 


— 8.962658 — 6.280293 


The effective mass is then recomputed from its exact ex- 
pression evaluating the integral by Simpson’s formula 
with eight intervals. The new value is found to be 
106.825m. The mode shape then compares better 
(white circles in Fig. 2) with the exact one and the fre- 
quency parameter is overcorrected 

162.558 

ut = 1808.74 = 2752.4 
106.825 


but lies within 6.26 per cent of the exact value. 


(9) CONCLUSIONS 


It is suggested that normal modes of complete air 
craft structures in free flight be computed on the basis of 
lumped inertia models. The results would appear to be 
reliable provided (1), the models be obtained as a result 
of some formula for numerical integration combined with 
collocation applied to the integral equation of the 
problem; (2), the constraints used to evaluate Gj,,. con- 
form with those of the static rig for comparison against 
experimental data; and (3), the higher frequency modes 
be corrected for effective mass and frequency by tech- 


niques similar to the one described. 
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y[x + (2/2)] (1 2), y(v) + v(x + h) + 


} 


1 +) [y’(x yi(x + h)]} 
The slopes may be computed in the isostatic reference 
frame since their differences only are involved. From 
ordinary beam theory we have for the slopes at the sta- 
tions due to the concentrated inertia loads 
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lz O + 2.484788 
a © Mx) = + 7.150795 
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TABLE 2 
res Numerical Results for First Elastic Mode from Lumped Mass Model 
he 
~ if 
he 
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nt — ‘ 5 =e |fo II 16 Ir 0 
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-65 I00 -I0 -20 -5 -~1369.394 
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A= Too -26 -20 I04 -44 -14 Xi = -~I540. 722 
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| TABLE 3 
| Numerical Results for Second Elastic Mod Siseas Rian Mass Model and Deflation of the First Mode 





-451.729 -339.348 -229.858 -619.324 -215.741 
1445.119 1177.520 639.779 1727-412 690.173 
1625.92I 1023.55I I0I19.657 2530.348 776.522 
-88.804 -144.588 49.174 338.628 -42.412 
-2804.629 -I767.308 -I806.955 -5273.646 -1339.460 


m a 
AG; soMA, = TIOS9D0ET 


x)= ( 1.170 186 -13.004 797 5.110 216 1.281 949 -25.951 253 ) 


Mi = I 608.74 x!M x,= 162.558 m. 
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APPENDIX 


Symmetrical Extended Coefficients as a Limiting Case of 


Dynamic Flexibilities 
Harmonic excitation p = ey produces a forced re 


twl 
> 


sponse gq = ex whose amplitude must satisfy (C 


w*M)x y. To the expansion of x in normal modes 
x 2. au 1 z, BX cr 
i 


36) 


there corresponds an expansion for y, that in view of 
Eqs. (4) and (5) may be written in the form 


tS Bla? —o)Mey, (37 
l 


n 
y w?) 8, Mu 

l 
So long as the forcing frequency is neither zero nor equal 
to one of the natural frequencies, the expansion (37), 


similar to Eq. (1S) is capable of representing any forcing 


amplitude. The coefficients must then have the values 
/ , 
Ui) V Vor) y 
a 5 , 6, = > , / 
wt, ;)" Mit (w,? — w*) Xi)’ Mx 


which, when substituted in Eq. (36), yield an expansion 
formula for the dynamic flexibilities matrix 
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For restrained systems the influence coefficients matrix 
is obtained as a limiting case by letting the forcing fre 


quency vanish 


It may be shown that the same expression applies to the 
matrix G for unrestrained systems. In other words, 


this matrix obtains as a limiting case of dynamic 
flexibilities by removing the free mode components of 
the expansion which have a simple pole at w? 0, be 
fore passing to the limit. 

To prove it, denote by G the expansion (3S) and ob 
serve that it already is symmetrical from structure and 
that it satisfies the relation 


AG G 


as a result of Eqs. (13). There remains to show that it 


satisfies Eq. (25). To this effeet multiply to the left by 


C and use Eq. (4) 


Now, on account of the expansion (20), there follows 


me Muti’ 
ise. 2S. «a 
1 Uy’ Mu 


Q.E.D. 


A Critique of Shock-Expansion Theory 


(Continued from page 680) 


{ ( WwW \/ ) I = :. 
Ww ; & Muy, 
* Li Xr’ 
Zz er 
; te” = @ Wx 
Mach Number of the main stream exceeds 2 7/3, and 


for Mach Numbers below this value the errors are 
additive. The fact that there is no cancellation of errors 
for low supersonic Mach Numbers suggests that at val- 
ues slightly above the critical value of 2/7/35 the process 
of balancing out of errors will be far from complete. 
However, this does not necessarily imply that the 
total error is larger for small Mach Numbers, since, 
in this case, the two individual errors are each much 
smaller for a given thickness ratio than in the hyper- 
sonic case. There is a possibility however that, for 
practical values of the thickness ratio, there will re- 
main a range of free-stream Mach Numbers of the 
order of three in which the individual errors have be- 
come large and the cancellation process has become 
insufficiently effective. Apart from this slight doubt, it 
would appear that shock-expansion theory is valid 
over the complete range of conditions likely to be en- 
countered in supersonic aerodynamics as distinct from 


“‘super-aerodynamics.’ 
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Reversal of Skin-Friction and Heat-Transfer 
Trends in High-Speed Laminar Flows 


jseorge M. Low 
Lewis Flight Propulsion Laboratory, NACA, Cleveland, Ohio 
April 6, 1955 


* NERALLY IT IS ASSUMED that, for all Mach Numbers, skin 
friction is increased by favorable pressure gradients, and 
decreased by adverse ones (see, for example, references 1, 2, and 
3 The purpose of this note is to show that, at sufficiently 
high Mach Numbers, the opposite is true—skin friction (and 
heat transfer) are lower on surfaces with favorable gradients 
than on a flat plate, and conversely 

The contentions of references 1, 2 and 3 are based on Falkner 
Skan types of solutions of the boundary-layer equations—i.e 
, = ox”. Results of these solutions are expressed in terms of 
the variation of a friction parameter as a function of a pressure 
gradient parameter (see Fig. 1, reference 2) However, the 
pressure gradient is contained implicitly in the skin-friction 
parameter; thus it is difficult to infer how the actual friction drag 
behaves as a function of pressure gradient 

A simpler relation between skin friction (or heat transfer) and 
pressure gradient is obtained from a perturbation analysis pre 
sented in reference 4. Admittedly, this analysis applies generally 
only when the pressure gradient is small. But the method can 
be applied for all pressure gradients (large or small) to obtain the 
rates of change of boundary-layer characteristics near a leading 
( Ilge 

The analysis leading to the present results is the same as that 
of reference 4. The basic assumptions are viscosity proportional 
to temperature, isothermal surfaces, and a Prandtl Number of 
1.72. It is further assumed that all the external flow character 


istics and boundary-layer variables can be expanded in Taylor 


series about their values at x = 0. Because emphasis is placed 
on rates of change of quantities near x = 0, only the linear terms 
in these series remain. Thus, the pressure gradient is character 


ized by a parameter Ps 
L f{du, 
P 
uy dx r=(0) 
where L is any length, wu is the velocity in the x-direction, 7 refers 
to conditions at x = (), and e refers to external conditions for all 


A drag parameter D and a heat-transfer parameter Q can then 


be obtained from reference 4, Eqs. (58) and (63 


| d , = (() } 
D -—— --—M, 2 
P dx r r=0 f"(0) 2 
1 fd q 1) 4,2, '(0) 
v = —- = UV? (3 


K s‘(0) 2 


P | dx Grp 1x =0 


lhe subscript fp refers to conditions on an equivalent flat plate, 
M, is the Mach Number at the leading edge, the functions g, f, H 


and s are tabulated in reference 4, and K is defined by 


ss l I t 
K- (‘«) , 
s(Q) t t ; 


(t,, and ¢,,, are the actual and adiabatic wall temperatures, re 
spectively.) 

Although both 7 and 7;, (also g and q;,) are infinite at x = 0, 
the ratio r/7;, (or g/gy») is one there. The rate of change of this 
ratio may therefore be regarded as a measure of the relative 
shearing stress a short distance from a leading edg¢ When D 
(or Y) is equal to zero, the rate of change of skin friction (or heat 
transfer) on an airfoil is equal to that on a flat plat When 
D is greater than zero and / is positive, the skin friction on an 
airfoil is greater than on a flat plate at the same short distance 
from a leading edge. When D is negative and ? positive, the 
converse is true 

The remainder of this discussion will be limited to favorable 
pressure gradients (P positive); the analysis and curves to be 
presented apply equally well to adverse pressure gradients. It 
is well known that, when J = 0, the drag and heat transfer on 
an airfoil are greater than on a flat plate; that is, for 1f = 0, both 
D and Q are greater than zero. However, there exist conditions 
of wall temperature and Mach Number such that the opposit« 
is true. The dividing conditions, that is the conditions for 
which the drag (or heat transfer) on the airfoil initially decrease 
at the same rate as on a flat plate, are given by Eqs. (2) and (3 
for D (or Q) equal to zero. These dividing conditions will be 
called the conditions for skin-friction or heat-transfer reversal 

In Fig. 1 the conditions for skin-friction reversal are shown 
To the left of the curve, the friction drag on an airfoil (P? positive 
is higher than on a flat plate; to the right, the drag is lower on 
the airfoil than on the plate. Even with maximum cooling (f 
0), skin friction on the airfoil will be greater than on the plate for 
VU < 1.67. For WV An airfoil at the 
temperature of an insulated plate (¢,,/,,,, = 1) has a drag equal 
to that of the flat plate at JJ = 4.71. Skin-friction reversal 


does not occur for temperature ratios greater than 1.47 


1.67, reversal is possible 


are shown in 


The conditions for heat-transfer reversal (VQ = 0 
Fig. 2. Heat transfer on an airfoil is greater than on a flat plate 
to the left of each of the curves shown, and conversely Note 
that, for cooled surfaces, heat-transfer reversal takes place at 
Mach Numbers of the order of one (The heat-transfer curve 
is discontinuous because only the heat transfer on a flat plate is 


zero when the temperature ratio equals onc On an airfoil the 
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Fic. 2. Conditions for heat-transfer reversal (Q = 0). 


rate of heat transfer cannot vanish when the surface is isother- 
mal.) 


Because the skin friction of the airfoil is referred to an equiva 
lent flat plate, the reversal trends are affected only by the ve 
locity and temperature profile shapes, and not by the surface 
values of viscosity and conductivity. A physical explanation 
for the skin-friction reversal is as follows: For any given Mach 
Number, the density along an airfoil (with a favorable pressure 
gradient) decreases in a streamwise direction both near the 
surface and at the outer edge of the boundary layer. From the 
equation of state it is evident that the density near an isothermal 
surface decreases faster than the stream density. Layers near 
the wall become therefore relatively less dense than those at a 
distance from the surface. However, the pressure drop is every- 
where the same; the ‘‘lighter’’ layers near the wall are therefore 
accelerated more than the ‘‘heavier’’ layers away from the wall 
This increases the velocity gradient O0u/Oy for small values of y, 
and explains the low-speed behavior of skin friction on airfoils 
The effect is accentuated by heating the surface, because layers 
near the wall are then even less dense. Simultaneously the de- 
crease (with favorable pressure gradient) of density level of the 
entire flow field tends to thicken the boundary layer; this factor 
leads to decreasing values of Ou/dy. At high Mach Numbers 
the latter effect becomes the dominant one, so that skin friction 
on an airfoil is less than on a flat plate. A similar reasoning can 


be employed to explain the heat-transfer reversal. 


Although the foregoing discussion was limited to the case of 
favorable pressure gradients, the reversal curves apply also to 
flows with adverse gradients. In this case, the reversal is from 
lower friction at low Mach Numbers, to higher friction at high 
Mach Numbers 
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Some Series Developments in Unsteady 
Aerodynamics 


Peter F. Jordan 
Dynamics Section, The Glenn L 
April 11, 1955 


Aartin Company, Baltimore, Md 


| ppecoggei AERODYNAMICS in the sonic and supersonic rang: 
are usually described by means of tabulated functions 
(Bessel functions, Fresnel or Schwarz integrals In applica 
tions—say in flutter calculations—complicated combinations 
of such functions occur. The object of the present note is to 


draw attention to the fact that power series representations of 
the same aerodynamics are sometimes unexpectedly simp! 
and are often more convenient for both analytical and numeri 
treatment. Power series are particularly advantageous 
electronic calculating equipment is used 

The sample series below refer to two-dimensional flow (wit! 
one exception). Their simplicity springs from two sources 

(i) The binomial summation, which appears if the basi 
tegral, Eq. (2) below, is developed in series form, can be pet 
formed by means of the relation! 


n ' 
Rac (" = <me + a): 
} > (m+a+tb)! 


m=0 \™ 


in which neither a nor } need be integral 

(ii) Two separate series arise at first if Eq. (2) is expanded, 
but these obligingly combine to form a single series 

We derive first a series for the lift distribution p(x), starting 


from the basic integral 


P(x) = 
» ex » 
2y ; 2ky 
J | [w’(x — ££) tr wU(x — Edie sles if e) as 
Ve le 0 k 
Qk 
“ / 
al Ve w ae | I 2 
kK { 
where 
P(X, t) = pi * pl x ef" local lift 
w(x, t) = vw(x)e'”! downwash (continuous in [0, x 


In addition, xc is the chordwise coordinate, c the chord; x 


at L.E.;w = tvc/v = 2ik. y = 1/M =a/v; «=1-y7 
We are going to consider elementary chordwise modes z,(x) = 
r : ; 1 ; 

x and denote the lift function due to z,(x) by p,(x). First, how 


ever, we consider elementary downwash distributions w,(x) = 
x, denoting the lift function due to w,(x) by p,(x), so that 





p(x) = rpr-(x) + 2ikp(x) 
Inserting in Eq. (2), with w(x) = w(x), the series (see reference 
9)* 
~] , 
wa —iu 
; gh optus J 7. 
Ph &" ¢ é yugjdaég / ’ . 
J0 0 n aT 3 
_ 1/9 
¥/2 
» 
n—Z 
* Recently restated in reference 3 The more general function 
Pt 
>, pe. P f gA (1 E)mH ee ~8UE Ty (yutdé 
0 
which appears in the treatment of supersonic delta wings (see reference +) 1S 
equivalent to the series 
i? ; ' ! 
P ty ' win +r 
A, wu, D ; 5 iu : x 
2 trA+uti!l 
n=0 
u—p 
2 i 
) ) t 2m)'» p 








The 
1 
For 
require 








READE! 
we ob fter some transformation, making use of Eq. (1 
“a ; 
\ N D ] ; 
2 2D 21 
D (ie i) 
\ “ 
> } 1 — 2? ! _ ] 


which is remarkably simple in its applica 


[his series* for pr(a 
‘1 ce the factors D,(.\7) are tabulated, is valid for any real 
ly ] 

value l 


xt the lift function p,(x) we describe the mode 


In discussing n 


f oscillation by 


X S ( \ S ( x 
Lao boas 
if) ) 
Thus we now assume ry = QO, I, 2 - we find that the cases 
nd ry > 1 have to be distinguished 


n+ 1)D, — 2D } 


A striking feature of the series Eq 1) is that the factor of the 


first power of k is proportional t« v), and independent of r if 


Suppose that we are interested in the stability at low 
frequency Tf of i given mode \ As 
\7? __ ) 
Dy — D D 
- 1) 
see immediately that for any real mode 2(a 
kD 
Im p(x (M2 — 2)sla y + OCR 5 
ly 7 l 
i 
rhe required stability criterion is given by the integral 
*] 
| VW? — 2)s(x) + 3(0) |s(x) dx 
JV 


rhe oscillation described by v) is stable (at low frequency 


ifA>0O 
For actual flutter calculations generalized force coefficients ar¢ 
ind define them by the rela 


equired. We denote them by C 


10n 


where Ly exp (vt) is the total lift per unit span, ZL; exp (rt) is 
and so on. The familiar 


the first order moment about L.E., 
quarter-chord flutter coefficients are related to the first four of 


these C,,, as follows 


It is easily shown that Eq. (3) fulfills the reverse flow relation 


t Recently investigated in reference 5 for the case of parabolic camber 


The higher coeflicients ( ire requir 
camber, compare references 3, 5, and 


Integration of Eq. (4) yields 


n=1 ‘ 
) 
( S k } 
wher 
n+] 2 
n+m--+ i ; 2 
} » 2)9 i. ¢ 4 
I 2/xr)D } 1/2)F 1/(n + 1)JE 
Separating real part and imaginary part 
yy + 9 
( k? >> (—k?)’ E,' 4 
— Dy 4 7) 4 » 
uv = _ 
™ 2n + 
k> k?)’ I 
—— on 4 » i ] 
) 
( Do (HRB) erm?" F k >> (-k } 
0 a] 
with 
I I I 
l ] l 
F I _ E I F I 
2 Zn + | 2 on + 2 
lable 1 gives a short tabulation of the coefficients EF,’ and E,,”.** 
A particular advantage of the series Eqs. (6) is that they hold 
no temptation to tabulate flutter coefficients over a 2k? + 


\7? — 1) instead of over k 
rhe series (6) converge, in principle at least, fi 
kand > 1 Admittedly, convergence is slow if & is large; for 


V7 ~ 1 and large values n 


[ + § with 6 small, and the forces 


However, if @ is large , then J = 
coefficients for ./ 1 should be used in preferen 


l + § (see reference | 


to the ‘‘cor 


rect’’ coeflicients for \/ rherefore, 


the series (6) remain advantageous for numerical calculations 
throughout the range required for practical flutter calculations, 


even with desk calculating machines heir advantage is even 


more pronounced if electronic equipment is used. As an ex 


ample, consider the diagnonal energy term which arises with the 


camber mode z,(x) = x(1 — 4 
ee & 
( = pix) — pola v — x? )dx Cy — ¢ — ( + ¢ 
" @¥ O 
_ n 
S* (-ik)" F,, (M 
— (n + 2) (n+ 4 


As the relatively few coefficients F 7) can be stored in the ele« 


tronic memory, the electronic calculator will, for a given value k, 


find the sum of this series in a shorter time than it would requir« 


to read a single value from a lengthy table of Bessel functions, 


say 
Series for the case 1/ 1 (see reference 1) are here added for 
completeness: to Eq. (3) corresponds 


' 
r'w,(x nal — ikx n— 05 1 + 2n 
1 — S = 7 
npenar n! n+r—0.5)!1— 2n 


VV 2rkx n =0 


en in reference 2 


rhe series for Co, Co, Cw and Cy; have already been giv 
Martin 


The Glenn 


** More complete tables can be obtained from 


Company 
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TABLE | 
M n=0 1 2 3 4 5 6 7 
1.25 E,,’ 3.0180 18.009 {8.185 742 72.804 50. 124 25.415 9.881 
EB,” 1.6077 8.3835 31.912 356 76.970 62.952 36.980 16.348 
1.5 EB,’ ©. 9111 2.1865 2.1568 1977 0.4274 0.1062 0.0194 0.0027 
E,,” 1.1388 1.6399 2.3615 7097 0.7518 0.2220 0.0470 0.0075 
S E,,’ 0.2450 0.3086 0.1428 0359 0.0058 0.0006 0.0001 0 
E,,” 0.7351 0.3267 0.2299 0764 0.0151 0.0020 0.0002 0 
lo Eq. (6) corresponds be a linear function of the absolute temperature 7. The analysis 


1-1 (—1k)" 2n — 1 2 
= 4 p » 19 9 ‘ 5 
rV 2ak 4 (n — 1)! 2n + 2m — 1 n— 1.5 
l1—1 “. (— tk)" 2n? + 1.5 ; 
( _ y! (8) 
e . ' 9 9 on — 
rV/ 2ak rr n! 2n + 2r + 2m l 
(n — 2.5)! 
(n+ r— 1.5)! 


The series Eq. (7) and Eq. (8) converge rapidly for all values k 
occurring in flutter calculations 

Finally, another example of the usefulness of such power series 
in analytical considerations: the limit at JJ = 1 — 0 of the 


Possio equation can be written: 


. 
x 
wx) = k | P(x — E\K(RENdE (9) 
J0 
with the kernel 
; 1+ (—ix)” 
K(x) = Ta ae 
oy 0 (n — 0.5) 


This power series representation is again remarkably simple 
(compare the representation of A(x) by means of Fresnel integrals 
as given in reference 7). By elementary algebraical manipula- 
tion it can be verified that pr(x), Eq. (7), is in fact the solution of 
Eq. (9).! 
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Simple Shear Flow Past a Flat Plate in a 
Compressible Viscous Fluid 


Ting Yi Li 
Department of Aeronautical Engineering, Rensselaer Polytechnic 
Institute, Troy, New York 


April 19, 1955 


N A PREVIOUS NOTE, the present writer discussed the simple 

shear flow past a flat plate in an incompressible fluid of small 
viscosity.!_ The purpose of the present note is to treat the same 
problem in a compressible viscous fluid. It will be shown that 
this problem in the compressible flow regime can be reduced to 
the corresponding problem in the incompressible flow case pro- 
vided that (1) the Prandtl number @ of the fluid is assumed to be 


unity and (2) the viscosity coefficient yu of the fluid is assumed to 


proceeds, in the conventional notations, as follows: 
Continuity equation 
(0/Ox) (pu) + (0/Oy) (pv) = 0 j 
Momentum equations 
pu(ou/Ox) + pi(Ou/Oy) = (0/Oy) |ulOu/oy 2 
(Op/oOy) = 0 3 
Energy equation 
pu(OlT/Ox) + pr(OH/ dy (0/Oy) [u(OH/oy (4 
where 
H = C,T + (1/2)u? 5 
represents the specific total energy of the flow In the derivation 


of Eq. (4), it has been assumed that ¢ = 1. The proper boundary 


conditions are 


©, HU = mw, (Ou/OY) = ta, T= 1 7 


where the subscript w denotes conditions on the surface of the 
flat plate and the subscript 0 denotes conditions pertaining to 
the undisturbed free stream evaluated at y = 0. The undis 


turbed free stream consists of a simple shear flow as follows 
Un = Uy + wy v. = 0 8)* 


where the subscript denotes conditions in the free stream 


In the present problem, uw, w), 7), and 7, are all assumed con 


stants. It can easily be shown that the solution for // is 
(H — H,,)/(Ho — H,) = u/u 4 
where 
BH, = €.7, Hy = (1/2)uo? + C,7 10) 


The continuity equation implies the existence of a stream fun 


tion ¥ such that 


(Oy Oy) = (p py il —(doy Ox) = (p/ po) 1] 
lf (p/p) = ( 7/7), the momentum equations lead to the follow 
ing equation for y: 

(Oy /O¢) (0?y/Ox0¢) — (O~/Ox) (O*Y/OF?) = wm (08/0 12 


where » = (uo/p)). In Eq. (12) the new independent variable ¢ 


is defined as follows 


ey 

5; = | (p/ po )dy 5) 
J0 : 

The boundary conditions become 
¢=0, (O~/ox) = (dy/oF) = O 14 
t— ©, (d¥/d¢) = mw, (d%y/de2) = +a 15 

Let 

n = (OV W/V wx) E = [wo/(uo/V m%X/ uo) 16 


and expand y in ascending power series of & as follows 


Y = V vwottox [foln) = Efi(n) +... Li 


* It can be shown that in the undisturbed free stream 1/2)u T 
Cpl H., along the stream lines It is assumed here that uo ( y=0 
/ (T o)y=0, Ho = (Ho)y=0 
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4s in reference 1, the ‘‘vorticity’’ number is — <1. Substituting 
Eq. (17) into Eq. (12) and collecting like power terms in — vield 
fofo” + 2f,’"’ =0 (18) 
Df 4 fee — fie? 4. ewe a & 19) 
vhere the prime denotes differentiation with respect to 7. The 
joundary conditions are 
7=0, fo =fo’ = 0, ff =f,’ = 0 20 
no, fy =1, fi” =1 21 


The problem of a simple shear flow past a flat plate in a com 
pressible viscous fluid is therefore reduced to the corresponding 
problem in the incompressible flow case which was treated in 
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Lifting Pressure on Delta Wings with Subsonic 
Leading Edges, Symmetrical Plan Form, and 
Discontinuous Slope 


Paul G. Kafka 
Mechanical Engineer, Boeing Airplane Company, Seattle, Wash. 
March 14, 1955 


SYMBOLS 


AP net lifting pressure 

free-stream dynamic pressure 1/2) pol 

free-stream velocity 
V free-stream Mach Number 

streamwise coordinate (parallel to } 

spanwise coordinate (normal to x 

angle between free-stream direction and line through wing 
vertex 


angle between leading edge and stream direction 


V Mo? 1 
BIgd 
Bl y/x 


wing slope in the y const. plane 
I 


DISCUSSION 

— LIFTING PRESSURE on delta wings whose slope is given 
by 

d Tx” | y l 

Lomax and M. A 


s been investigated by H Heaslet in NACA 


IN 2497 They obtained the following result 
ae ee ; 
Ap x\? ra q be 
q B pT were 


Coefficients 6; are found by equating coefficients of like powers 


of » in the following expression 


iii , (n — tj)? tay 1— ? 
1p + q)! J1 


rhe symbol f / means finite part of the integral. 

The Lomax-Heaslet theory does not cover the important cases 
when p # (0 and q is an odd number, such as for instance \ = 
Ixy. It is the purpose of this paper to extend the theory to 
such cases 


As in the Lomax-Heaslet theory, the fundamental integral 


equation 


FORUM 725 
Mx, y) = 
Ap 
7 . x x x,y 
| 
/ d y ax 1) 
tr Jy re (7 — AnrY we — 4 — p(y — ¥ 


is inverted under the assumption that 
Ap/q C(@)y q fa 
For symmetrical plan forms this leads to the equation 


g\r ta tkeotia 
A= - 6° * 2 C(@)H(0, n)\dé 5 
x tr a 


« 


wher« 
y7— il ‘ l1—f? 
| \ dt —#, 80 , 
6- ft ate, 
H0,9) = \y8, > 
(n —? in/i—? 
at a a 
J-1 (0-12)? re 
& = Btgsd 
Eq. (5) is an integral equation for the unknown function C(@ 


If \ has the form (1), Eq. (5) can be reduced to the form 


I1(6, n)dé 7 


l derivative ol n*, for odd >. inishes ever\ whe re 


The p + q+ 


except at » = 0, where it behaves like a multipole singularity 
Ss 44(n 
rp! p+ 
2p — | Re b+ | 
S n) = lim —- % 
al r(n? + u i , 
=] 2 uy S 
r 1,3, 5 
Thus Eq. (7) can be written as 
l Yi — 70" nm gr*e C(O 
Taf a (7) = — dé g 
. tr n On” *4 F 9 ” 


00 gi aC(@ 
If we put dé = y(n), Eq. (9) isa differential equation 
= @—7 


for ¥(n It states that the p + g + 1 derivative of ¥y must have 
a singularity not stronger than 7S, n)atn (), and must vanish 
forn +0 


For positive n, the most general solution of Eq. (9) which satis- 


fies these conditions, is 


74, + . + 
erp *aC(@ i am 
| dé = > a et S acn' - 10 
a a a 
uv 6 6 7 0 q+ 
where 1(m) is the unit step function. Inverting integral Eq 
(10) we obtain 
I ity V Ot — 7 5 Y 
er +a C(0) = { — — le — 
T-V O° — O e 6 


+] F dn 11 


where A is an arbitrary constant Evaluating the integrals this 


gives 
p+aqt+l b.e* Prq " 
er *aC(@) = ) ) C0 Ch 12 
» 9 ° @ 
0 ye" - & q+l 


and finally, introducing Eq. (12) into Eq. (4a), we obtain the 


lifting pressure 





Ap x aj Ptatl 50 cee 0 
( ) + + doy 0° Ch las 
q 8 ee: V A? — 6 71 of 
0, 2,4 


where 7 has been restricted to even values, since the pressure must 
be an even function of 6 
Coefficients ); and c; depend only on 6, and will be determined 


presently. Substituting Eq. (12) into Eq. (5) we obtain 


q+1 


p+qtl1 soe Bae) b0 
Tn 9 = > 
47 3 9 ri) V 072 — & 
Pr ¢ A, °7 | — ;)pta i — 7 
> 6:0! Ch ) ao f | Li V dt 
q+ Ji (@ — tyr tat? 
n a) 14 


After some lengthy calculations, this leads to 


J! f 
p+ q+ l biti p+a 6 
p> = . f Z. ot! sin7! ; dt, n 0 
0 V tt — 6 oat 
1=0, 2. } (15 


The 0;'s and c¢;’s are found by equating coefficients of like powers of 


n on both sides of Eq. (15 
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Frictional Heating of Nonisothermal Walls 


Myron Tribus and James E. MahImeister 

Engineering Department, University of California 
Los Angeles, Calif 

April 19, 1955 


IT A LETTER TO THE READER'S ForuM, Korobkin has suggested 
an experimental method for determining the “effective wall 
temperature”, 7.s, and the “‘heat-transfer coefficient’, . We 
would like to give our interpretation of the significance of the 
measurements made in the manner suggested in his letter 

In reference 2, it was shown that the heat transfer from a 
body subject to aerodynamic heating is given by Eq. 1 


f=x — va 
q(x) - ff ‘. h(é, x)d[7.(&) — Taa.(é)] (1) 


Eq. 1 is to be interpreted as a Stieltjes’ Integral. g(x) repre 
sents the heat flux per unit area at a point x on the body. 77, 
represents the wall temperature; 7,4. represents the adiabatic 
wall temperature. The symbol — is a dummy variable which 
takes into account the effect of the upstream condition on the 
temperature at a point x. A(g, x) is an integrating kernel 
described in reference 2 and also in reference 3. The heat trans 
fer at a given point can be made to be equal to zero either by 
having 7,, equal 7,4. everywhere over the surface of the body 
(true adiabatic wall temperature everywhere) or by so adjusting 
the wall temperature that the local heat transfer is zero. In 
reference 1, it is suggested that the wall temperature be made uni- 
form and so varied to produce zero heat transfer at a particular 
point. Under these circumstances, Eq. 1 may be integrated to 
give for uniform surface temperature Eq. 2. The symbol 7ya.(() 
refers to 


- dT 4a.(é) 


dé 


g(x) = [T,, — Taa.(0)]h(0, x) - h(é, x) dé (2 


u 
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the adiabatic wall temperature at the stagnation point. [t 
therefore equal to the stagnation temperature The unifor 
wall temperature required to make equal to zero is therefor 


given by Eq. 3 


ro I dT 4.(3 
lott. = 1. = Tei. 4 h(t, x di 
h(O, x 0 di 


Eq. 2 reveals that the heat transfer should vary linearly with t} 
wall temperature, and that the slope of the heat flux versus w 
temperature line is the ‘‘heat-transfer coefficient’’ (0, x) appro 
priate to an isothermal surface Eq. 3, however, shows that 
the ‘effective wall temperature” obtained from the graph of 
versus 7), does not necessarily go through zero at the local 
abatic wall temperature. If in the practical application of the 
data gathered by the method of Korobkin one is dealing with 
isothermal surface, it can be seen from Eqs. 2 and 38 that th 
combination of the “effective temperature’ and “‘heat-transfer 
coefficient’’ obtained by his method will yield the correct answers 
On the other hand, if the surface is to be nonuniform in temper 
ature, both the effective wall temperature and resistance to heat 
transfer will be in error 

Data gathered by the method suggested by Korobkin can by 
used to find the true adiabatic wall temperature provided one or 
another of the integrating kernels of reference 3 applies 
Tut. is the effective temperature as defined by Eckert* and this 








temperature, along with /(0, x) is determined by experimen 


then one has, from Eq. 3, 
h(O, x) [Teors.(x) — Ts 


Using the methods of reference 2 one may invert this integral 


equation and obtain 


Yeu, ~ Tat * { "g(t, x) W(O, 2)[Tataz. — Tett.(t)]dé 
¢ 0 


The kernel-pairs ¢(£, x) and A(é, x) are given in reference 3. Eq 


(5) is easily handled graphically 
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Direct Calculation of Propeller Deflection 


S. M. Ramachandra 

Department of Aeronautics, Indian Institute of Science, Bangalore 
South India 

April 20,1955 


SUMMARY 


A simple noniterative method for the calculation of propeller deflection | 
obtained by a finite difference approximation of the differential equat 


governing the bending of the propeller 


NOTATION 


CPF; centrifugal force at the point x; 
E; Young’s Modulus of elasticity at the point x; 
h the spacing of the points x;, assumed constant 








obtained 
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moment of inertia of cross-section at x;, about a line passing ~ 
It through the center of gravity of the section and parallel to the r i - 
lifor; chord , j 4 pa 
refor vi the bending moment at the point a j i _ te “ , 
thrust force at the point x; ee | ee —TrE p gS ee ¥ 
distance of the 7th point from the root of the blade 6—— fe bs 
f t +1 
leflection of the point x; ee ee ee 
INTRODUCTION ’ _ 
, q HE PROPELLER BLADE 1s subject to the action of the aero- y - 
i) > ° 
dynamic forces and the centrifugal forces Due to the aero , - 
' forces the blade bends A first-order estimate of the 
t] lade deflection can be made by ignoring the variation of the 
1 Lae : ‘ : oa IK 
| verodvnamic forces with the deformation of the blade The 
oO 
lade iv be treated as a cantilever subject to axial and trans 
ay verse loads Assuming small deflections, the bending equation 
Pi ; P WEI , With this in view the blade is divided into n equal parts, each 
h a-y/dax 1 4 ) : = 
or of width equal to 4. By a finite difference approximation 
tl ° » > 
sf may be used. Because of the complex nature of the forces and 
ster : d*y/dx? 1 /n* y — 24 y 2 
; the flexural rigidity, an approximate solution of this equation is ; . } 
i obtained by determining the deflections at a finite number of at the 7th point The bending moment at the point x; (Fig. | 
1 
2 points is equal to 
1€ 
1 be 
Ce O 
mr WV, = h{7 + 27 + n—1)7 — CF 4 - yj Cr ) - CF,(4 } 
it The Eq. (1) then becomes 
wT. as 7 he [ 
24: + ¥ + CF ¥ — } + CF 1 —y¥ CF — y ; ] 27 + 
Ej! El 
, 1 
ra 
which may be rewritten as 
- h? ad ; —_ } ; 
; 712+ — CF + CF + CF, ty 1 + CF + 
Ei] E;l 
( | j a h 
— CF a + + — Cry 1 + 27 + (n —1)7 4 
E;l E xl E;1 
rhis equation may be written for the first » points, by giving 7 the values L123 n Since vo, the deflection at the root 
is zero, the following n simultaneous equations are obtained 
h? te he h 
v{ 1 4 CF) + ye © CF; 4 + y, — CF, = — (1 + 272 4 + wT.) — 
. Eyl Eyl Eyl Ful 
“ h? h? \ h? h? 
uj2+ — CF, + + CF ] . lL. —— CF ) + 3— CF. y CF 
FE, l | E,l E,l E,\l 
A 
T, + 27 (a —1)8 
El 
h _ : he ee h? h ? 
— Kis +. -- CF; + CF + y ba =~ £2 + Ws CF; + + 9, =— OF IT; + 27, 4+ ) 
L E.] | E.] Fol E.] EI] 
+ + (n = 2)7 
h* a h? h 
= % 2+ = CF + y | + — Ch = — 1 
a | | a | | en | 
which may be written concisely in the form of a matrix equation 
iY=B8B 6 
where the column matrix 
‘ ) op 
Y= iy ¥ Vut r 
is the deflection matrix. (The matrix notation of reference | is used here 
' 97; + 272 + + nT, T2 + 273 4 + (n—1)7 1 | 
B h : , (8) 
| Eol El I I { 


ind the matrix A is given by 
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(6) is seen to be 


y = A~'B 


The solution of the Eq 


where A ~! is the reciprocal matrix of the matrix A, which car 


obtained by the standard methods given in reference | 
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Approximate Solutions of the Incompressible 
Laminar Boundary-Layer Equations for a Flat 
Plate in a Shear Flow 
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Assistant Professor of Aeronautical Engineering, Rensselaer 
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April 23, 1955 


SUMMARY 


The two-dimensional steady boundary-layer problem for a flat plate in a 
shear flow of incompressible fluid is considered. Solutions for the boundary 
layer thickness, skin friction, and the velocity distribution in the boundary 
layer are obtained by the Karman-Pohlhausen technique. Comparison with 
the boundary layer of a uniform flow has also been made to show the effect of 


vorticity 


SYMBOLS 


a4, ¥ Cartesian coordinates with origin at leading edge of the plate and 
x extending in the direction of the free-stream velocity 

%, 2 velocity components parallel to x and y axes, respectively 

l Uo + ay, velocity distribution of the shear flow 

Uo velocity of the uniform flow 

Sz pl'e/pai?x 

p density 

“ viscosity 

a vorticity of the shear flow 

T skin friction of the shear flow 

T. (7/pU vai*x), skin-friction coefficient 


ro = skin friction of the uniform flow L’o 
ot) boundary-layer thickness of the shear flow L 
boundary-layer thickness of the uniform flow l 


59 = 
INTRODUCTION 


igo BOUNDARY-LAYER problem of a flat plate situated in a shear 
flow of an incompressible fluid has been treated in referenc« 
1 by the method of asymptotic expansion, and some results per 
taining to the characteristics of the boundary layer of a shear 
flow have been presented. In this paper, approximate solutions 
for this problem are obtained by the Karman-Pohlhausen tech 
nique.’ It is seen that the velocity profiles at different stations 
on the plate are not similar to each other, and a form factor which 
depends on the boundary-layer thickness is needed to define them 
Since no single similarity parameter can be used to our advantage 
here, the boundary-layer thickness and skin-friction distribution 
along the plate are shown in graphs in dimensionless parameters 
from which the velocity distribution can be easily calculated 
For comparison, the boundary-layer thickness ratio and_ the 
skin-friction ratio between a uniform flow and a shear flow have 


been obtained and shown in graphs 


Basic EQUATIONS 
The two-dimensional, incompressible, and steady boundar) 
layer over a flat plate situated in a shear flow defined by 
U = Uy + ay I 
In Eq. (1), Uy and @ are constants and the un! 
form velocity Uy is always positive, but @ may be positive 0 


is considered. 


* The author wishes to express his gratitude to Vural Oskay for his help 


in carrying out the calculations 
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iG Variation of a)é6 and r. with S oO positive 


negative Phe boundary-layer equations without pressure 


gradient are* 


u{Ou/OxX) + v(Ou/OY) = v(O07u/OyV* 2 

(Ou / OX (O7v/Ov) = 0 =) 

where v = p/p Let a boundary-laver thickness 6(x) be defined 
such that uw = U' at y = (4 The boundary conditions are 

Q=v=0 O7u /Ov* = 0 aty = 0 j 

Mu l Ou/OV = a O7u /Ov? = UO ity = d(x 5 


By integrating Eq. (1) across the boundary-laver thickness and 
making use of Eq. (3) and appropriate boundary conditions, the 


following inomentum equation is obtained 


%O %o 
, Ou Ou- T 
l - a0 = ay = va h 
Jo OX Jo OX p 
where 7 u(Ou/Oy) at y = 0. The velocity profile # in the 
boundary laver can be expressed as 
“= l + a(x )| f(n ‘ 


with » vA Following Pohlhausen’s idea® let us take 


(n) = an + by? + cn*® + dn s 


If the velocity profiles are similar, the function /()) will have coi 


stant coefficients, a, b, c, and d. However, from the boundary 


conditions (5), it can be easily seen that this is not the case, aud 


i(n Fi n ab l + «wb)\G(n y 

where Fin 2n 2n*® + n 10 

and G(n n — 2n*? + 9 11 

In Eq. (9), ad “(0 + a6) is then the form factor for the velocity 


profile 


It can be seen that, when () is zero” ie., only simple shear 


flow is present f(y) reduces to yn and uw = ay. In this case, the 
solution of the boundary-layer problem is simply f(7) = n, or 
(= ay. This is due to the fact that the simple shear flow satis 


fies the boundary-laver differential equations as well as all the 
boundary conditions. Therefore, it is expected that F(n) — G(n 
should always be equal to 7, when a proper form for f(7) like Eq 
S)ischosen. If an additional term en is included in Eq. (8), it is 


found that 
Fin 5/2)n — An* + Jn 3/2)n 
(s(n = (3/2)y jn + An (3/2)n 
Hence, Eq. (7) ean be written as 
u UV F(n) + and(x = (\F(n) + ay 7a 
It follows that the skin friction 7 is given by 
r= pa + wl FO) /6(x 12 


Therefore, in Eq. (6), (7/p) — ve is simply equal to vl) F'(O) 


FORUM 720 
SOLUTIONS OF THE BOUNDARY-LAYER PROBLEM 


With the 


easily integrated to give 


velocity protile given by I ] v Kq 6) can be 


0.005 5556 O O24 36506 A 133 


In obtaining the above expression, the boundary condition (4 


(hata 0 has been used Eq. (13) can be written a 
O05 SEA a6 (6.024, 365( a6 pay*x  p 1/S 14 
vhere S; is a non-dimensional parameter and a a ly) has the 


dimension of L In reference 1, a similarity parameter for the 


/ 


boundary-laver equations defined as v6 Ly /Bvaqx was 


used. Eq. (14) shows that this is possible only for a)é sufficiently 
large 
A non-dimensional skin-friction coetticient 7, will be used and 


is defined as 


Phe skin friction 7 can be obtained from the expression 
T rT. p LU way*n 16 


Eq. (14) can be solved for a6, but it seems to be useful to plot 
i graph of a6 vs. S, for values of a)6 and S; which are com 
This graph is shown in Fig. | and can be used to 


In the 


monly used 
calculate the boundary-layer thickness for any given flow 
is also given, from which the skin friction 
For S; > 107, both 


curves are practically straight lines in the logarithm scale with 


same figure, rt, vs. .S; 
can be calculated by making use of Eq. (16 
slope equal to —(1/2 

For a given aq, the velocity profile « at any x can be easily de 
termined by making use of the curve a6 vs. S; and Eq. (7 No 
difficulty will arise as long as a; is positive and 6 is defined. If a 
is negative, however, one complication will arise. As soon as the 
value of a)6 becomes smaller than —3 /5, it can be easily shown that 


f(m) will exceed 1 at » = 1/2. Therefore, special attention is 


needed for that part of the plate where —a)é may exceed 3/5 
[his situation is similar to that of a boundary-laver problem with 
pressure gradient 

The boundary-layer thickness parameter a;é6 and the skin 


friction coefficient 7, for negative a; are plotted in Fig. 2. For the 
same S;, both —a)é and —7, for negative a; are smaller than 


than those for positive a;, but the difference becomes smaller and 
smaller as S; increases. They are practically equal to each other 


when S,; = 10) 


COMPARISON OF THE BOUNDARY-LAYER CHARACTERISTICS OF A 
SHEAR FLOW WITH THOSE OF A UNIFORM FLOW 


If a; is equal to zero, Eq. (13) gives the usual boundary-laver 


thickness for a uniform flow / 
b) = 5.835,6 V px/pl 17 


(he skin friction ts given by 
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J } for 4, posite 
&r ~ t Le + + + 4 
# C 6 | 
| 13 A \5 i> 
7\ r= . }/ 
vA 0G Sz : 


iG. 3. Variation of &/6) and ro/r with S, 


tT = 0.3438uV ply? /ux (18 


A comparison of the boundary-layer characteristics of a shear 
flow with those of a uniform flow can be made by examining for a 
given Reynolds Number pl’x/y the boundary-layer thickness and 
the skin-friction ratios at the same x. In view of Eq. (17), Eq 
14) can be transformed to give the thickness ratio as 


60/6 = 5.835,6(0.029,365 — 0.005,555a);6)' 2 (19 


The skin-friction ratio is given by 


= (20 
T 6, 2 + a6 

The calculated results are plotted in Fig. 3 for both cases, a; posi 
tive and negative. All the ratios will approach one as a; — 0 for 
any finite x. Of course, if a; is positive and becomes very large, 
6/69 may be so large such that the boundary-layer theory may 
case to be valid. In fact, Eq. (19) shows that (6)/6) = O when 
@6 = 5.286. Therefore, within the limitation of the present ap 
proximations, it seems that no valid solution can be expected for 
@6 = 4 


NUMERICAL EXAMPLE 


As a numerical example, let us consider the boundary layer of a 
flat plate with Reynolds Number = 10° and a, = 10 ft.~! at 
v = 2 ft 
From Fig. 1, the value of a6 is 0.117 


First, from the given data S, is found to be 2.5-10% 
Hence, the boundary- 
layer thickness at x = 2 ft. is 0.0117 ft. The skin-friction 
coefficient 7, = 0.0072 can be obtained from Fig. 1, and the skin 


friction 7 can then be calculated by Eq. (16). The values of 


6)/6 = 0.988 and ro/7r = 0.952 are obtained from Fig. 3 

It is of some interest to find the variation ot 6)/6 with x. If 
pUox/pu is also taken as 10° at x = 2 ft., the calculated results 
are given in Table 1. It can be seen that as x approaches 0, 
the ratio 69/6 will approach 1. On the other hand, as ¥ increases, 
the ratio will become smaller and smaller. Therefore, the 
boundary-layer thickness for a shear flow increases at a rate 
faster than that of a uniform flow with the same pl yx/yu, and 


hence at a rate greater than x /? 
CONCLUSIONS 
Solutions tor the boundary-layer problem of a flat plate situated 


in a shear flow of an incompressible fluid have been obtained by 
the Karman-Pohlhausen technique. Results are given so that the 


TABLE | 
pl ox/p 106 at x 2 ft 
x 0.01 0.1 0.5 1 2 4 S 12 
RN, 5-10% 5-104 2.5-108 5-10 10¢ 2-106 +. 10° 6. 10° 
Sz 5-1 5-104 104 5-16 2.5°10 1.25-10 667 411.75 
mn /b 0.999 0 996 O 994 0.991 0.988 0.983 0.977 0.900 


OCTOBER, 1955 


boundary-layer thickness, the skin friction, and the velocity dis 
tribution for « in the boundary layer can be found by simp} 
calculations. Comparison with the usual boundary laver of yy 
forin flow with a velocity l%) equal to that experienced by thy 
leading edge of the plate in a shear flow shows that for the sany 
Reynolds Number pl yx /u, the boundary-layer thickness and skj 
friction for a shear flow wih positive a; are greater than those fo 
a uniform flow When a is negative, the opposite is true. |} 
may be noted, however, that in both cases the difference is no; 
more than 8 per cent for S, > 10°. As the absolute value of q 
becomes very small (0.001, say) both ratios will approach on¢ 
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An Approximation to the Lift of a Two- 
Dimensional Cascade of Airfoils 


B. W. McCormick* 
The Pennsylvania State University, Ordnance Research 
Laboratory, University Park, Pa 


May 9, 1955 


I’ A SINGLE two-dimensional flat-plate airfoil is replaced by a 
single vortex at its quarter-chord point and the strength of the 
vortex adjusted so that the resulting flow is tangent to the airfoil 
at the three-quarter chord point, then the lift of the vortex will be 
equal to the lift of the airfoil as determined from the exact solu 
tion obtained by conformal mapping. This is the basis for Weiss 
inger’s approximate L-method solution of finite wings describe 
in reference 1. In order to examine the feasibility of extending 
the L-method to a cascade of finite lifting surfaces, the present 
investigation was undertaken whereby the validity of replacing 
each airfoil in a two-dimensional cascade by a free-vortex at its 
quarter-chord point and satisfying the aforementioned condition 


at the three-quarter chord point is established 


* Associate Professor of Engineering Research 
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FIG 2 
an airfoil in a cascade 


Consider the cascade of flat plate airfoils in Fig. 1. The origin 


of the coordinate axes, x and y, is taken to be the quarter-chord 


point of one of the airfoils 


FORUM 


Each airfoil is replaced at its quarter-chord point by a vortex of 
strength T. IT is now determined from the stipulation that, to 


the approximation of a being small, 


V=a 
where 
w = downwash induced at three-quarter chord pomit 
V = undisturbed velocity of free stream 
a = angle of attack of airfoil relative to | 


From the geometry of Fig. 1, 


2rwt ad > ncos 8 + c/2l ' 
¥ Ps n? + n(c/t) cos 8 + (c/2t 
where 
¢ = airfoil chord 
t = distance between airfoils as shown in Fig. | 
3 = angle as shown in Fig. 1 


his sum can be evaluated by the contour integral method of 
reference 2 
Using the relationship 
r= (1/2)cC, 1 3 


the section lift coefticient, C), can be obtained as 


( 2 t/c tan? rX¥ + tanh? 7} ‘ 
2 ra 7 sin 8 tanh r}(1 + tan? rX) + cos 8 tan 7X (1 — tanh? wr) 
cos 8 sin p 
where { = = 
2(t/¢ 2(t/c 


The quantity C)/2ra has been determined exactly by con 
formal mapping for this case and can be found in reference 3 
and 90°, the above result re 


For the two limiting cases of 6 = 0 


duces to 
C, tan } w/[2(t/c)]} 
= 8=0 
2ra w/|2(t/c) 
le 
C; tanh } r/[2(t/c)]} 
= 38 = 90 
2ra w/|2(t/¢ 


rhese results are identical with the exact solution. The agree 

ment for intermediate values of 8 between the exact solution and 
the approximate result (4) is presented in Fig. 2. The approxi 
mation is in close agreement with the exact solution over a wide 
range of values and approaches the exact solution for increasing 


values of t/¢ 
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Comment on ‘‘Wedge Pressure Coefficients in 
Transonic Flow by Hydraulic Analogy”’ 


H. J. Bomelburg 

Post Doctoral Fellow, Institute for Fluid Dynamics and Applied 
Mathematics, University of Maryland, College Park, Md 

May 11, 1955 


yy A RECENT PAPER Fleddermann and Stancil’ compare their 
own measurements of pressure distributions on wedge profiles 

WT 1L 4 ° +4 

1 tne transonic region obtained by the hydraulic analogy with 


corresponding measurements obtained by Griffith? in a shock tube 


- Guderley Mp1 
Jt ft 


\ Griffith Mpf 
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Pressure coefficients for a 2° double wedge for 2 
to this case by means 


Fic. 1 
Griffith's measurements are transferred 
of von Karman’s similarity law for transonic flows 
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and with theoretical calculations by Guderley.* They found out 
that the water analogy vields values that are by no means in a 
quantitative agreement with theoretical predictions or with the 
corresponding experiments in air. The author of the present 
comment had already made the same comparison in 1952, but he 
found a much better agreement.‘ Fig. | shows a diagram from 
the author’s paper that corresponds to Fig. 2 in the paper of 
Fleddermann and Stancil. How can this difference in the results 
be interpreted ? 

It is a known fact that the quantitative applicability of the 
water analogy is being appraised today by different investigators 
in a quite different manner, that ranges from a rather pessimistic 
to an Overoptimistic one. The author has shown in his cited paper 
to a greater leigth that it is by no means an easy experiment to 
run a shallow water channel in order to obtain reliable results, as 
there are so many possibilities which can cause errors in the results 
of the water analogy. For instance even the normal dust that is 
present nearly everywhere and that accumulates on the surface 
of the water affects the reliability of the water analogy more or 
less. Therefore the author has used kerosene (1/4-in. depth 
in his later experiments as the surface of kerosene remains always 
clean. Furthermore, it is just difficult to get reliable data from 
the water analogy in the transonic region owing to the existenc« 
of weak shock waves in this region as especially the analogous 
weak water jumps do not behave as single unsteady jumps but as 
series of smooth surface waves.® It is not possible to go further 


into details in this short comment that has only the purpose of 


[CAL SCIENCES OCTOBER, 


pointing to the more comprehensive papers on this matter. The 
iuthor might refer in this connection also to the extensive investi. 
gations of Crossley and Harleman® (cited already in reference 1), 
that vield substantially the same exact results apart from the 
well-known ‘‘overshooting’’ at the tip of the wedge that is due te 
the greater depth (1 in.) they used for their experiments 
Summarizing, it can be said that the results presented by 
Fleddermann and Stancil cannot stand as typical results from the 
water analogy as they ure presented. The explanations that are 


given by the authors for the discrepancies cannot be upheld 
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(tgation 


Design of Two-Dimensional Continuous-Curvature Supersonic Nozzles 


Continued from page 692 


CONCLUDING REMARKS 


The method of Riise and Puckett for obtaining con 
tinuous curvature in the design of supersonic nozzles 
is given herein. Also presented is a method for obtaining 
the length of the entire contour. Once the right com 
bination is found for the parameters 6, and A 64, the 
contour of the terminal curve can be obtained with 
automatic computing equipment. Such equipment is 
available if the value of #;p 0, is some multiple of | . or 
|» between the values of |» and 1. This corre 
sponds to a value of A 6, being some multiple of ' , or 
' « between the values of | and 4. This limitation of 
Op 0, =! 
as the rate of change of curvature becomes excessive 


,or A 6, = | isa practical one, inasmuch 


near the characteristic point as #;p 6, approaches | 
(A @, approaches zero), thereby greatly reducing the 
abilitv of the elastic curve to match the theoretical 
aerodynamic curve. Over the range of Mach Numbers 
for which the nozzle ts designed, it is desirable that both 


4, and A 4, be continuous functions of the Mach Num- 


ber. 
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